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1  Introduction 


This  report  describes  research  on  computational  aspects  of  control  and  stabilization  of  flexible 
structures. 

The  focus  of  this  work  was  on  control  observation  and  stabilization  of  elastic  beams.  Two 
models  were  explored,  the  Euler  Bernouli  beam  model,  and  the  Timoshenko  beam  model. 
In  both  cases,  the  computational  approximation  method  chosen  was  that  of  finite  elements, 
piece-wise  cubic  in  case  of  E-B  beam,  and  both  piece-wise  cubic  and  piece-wise  linear  in 
the  case  of  Timoshenko  beam  model.  The  immediate  objective  of  this  research  was  to  gain 
insight  into  the  dynamics  of  control/observation  process  through  numerical  computation,  as 
a  complementary  tool  to  analysis.  This  included  both  linear  output,  and  state  feedback, 
as  well  as  more  advanced  control  concepts  such  as  adaptive  observer  and  adaptive  model 
reference  control.  The  software  developed  in  this  work  will  serve  for  exploration  of  flexible 
structures,  composed  of  beams  and  some  problems  involving  nonlinear  phenomena  arising 
in  control  of  large  flexible  structures. 

Currently,  the  established  practical  methodology  for  simulation  and  control  of  large  flex¬ 
ible  structures  consists  of  developing  a  large  finite  element  model  directly  from  mechanical 
properties  of  the  structure  represented  by  a  finite  system  of  masses  and  springs,  such  as,  for 
example,  the  approach  underlying  the  NASTRAN  code.  In  this  work,  we  proceed  from  a 
continuum  model  of  the  structure,  described  by  partial  differential  equations,  to  the  finite 
element  model  as  an  approximation  of  the  continuum  model.  By  examining  a  family  of 
finite  element  models  of  increasing  dimension,  one  can  in\3stigate  the  control /observation 
phenomena  related  to  the  approximation  of  the  continuum  model  by  the  finite  dimensional 
finite-element  models.  This  approach  offers  a  possibility  of  gaining  insight  into  the  rela¬ 
tionship  between  continuum  and  finite  dimensional  models  of  flexible  structures,  and  into 
questions  such  as  the  role  of  neglected  high-frequency  components. 

Currently,  there  exists  a  large  gap  in  research  methodologies  between  two  approaches: 
the  PDE-based  theory  of  dynamics  of  elastic  structures,  and  the  more  practice-oriented 
methodology  based  on  finite  dimensional  models  derived  from  structural  mechanics  models. 
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The  corresponding  control  approaches  are  very  far  apart.  The  present  research  is  aimed  in 
part  at  providing  a  link  between  the  two  by  using  computational  tools. 

2  Discussion  of  Models 

Traditional  linear  PDE  models  of  elasticity  such  as  Euler- Bernoulli  (E-B)  equation,  and 
Timoshenko  equation,  are  derived  under  certain  assumptions  of  linearity  of  stress-strain 
relations,  and  assuming  small  amplitude  of  deflections.  These  models  have  an  infinite  number 
of  eigenvalues  placed  on  the  imaginary  axis.  Characteristic  frequencies  associated  with  these 
eigenvalues  from  an  increasing  sequence  with  nc  bound.  In  modified  E-B  models,  a  viscous 
damping  is  introduced  resulting  in  the  uniform  shift  of  the  eigenvalues  to  the  left  half  plane. 
This  property  is  inconsistent  with  the  experimentally  observed  fact  that,  in  flexible  structures 
the  modulus  of  frequency  characteristics  eventually  rolls  off  at  higher  frequencies.  This  fact 
has  led  to  several  efforts  in  constructing  models  that  provide  frequency-dependent  damping. 
Except  for  the  Rayleigh  damping,  these  are  either  nonlinear  models,  or  models  with  materials 
memory  governed  by  PDE’s  with  "hereditary”  terms.  Such  models  may  be  more  accurate  at 
higher  frequencies,  but  they  introduce  additional  complication  such  as  nonlinearities  or  an 
augmented  state  space.  It  has  not  yet  been  proved  that  the  control  design  problem  benefits 
from  introduction  of  these  additional  complications. 

In  this  report,  we  still  rely  on  linear  PDE  models,  being  aware  that  results  concerning 
high-frequency  behavior  of  such  models  may  be  not  representative  for  the  control  design. 

An  additional  loss  of  accuracy  at  the  high  frequency  range  is  introduced  by  the  finite 
element  approximations.  It  is  known  that  for  both  Euler-Bemoulli  equations  and  Timo¬ 
shenko  equation,  the  eigenvalues  of  the  approximate  finite- element  models  do  not  all  match 
true  eigenvalues  of  the  PDE  model.  In  both  cases,  only  the  first  few  modal  frequencies  are 
reproduced  with  negligible  numerical  error  (say  less  than  0.1%),  and  the  mismatch  increases 
for  higher  frequencies.  This  mismatch  is  illustrated  numerically  in  a  subsequent  section  of 
this  report.  Hence,  for  high  frequencies,  not  only  the  continuum  models  are  inaccurate,  but 
their  usual  computational  approximations  are  inaccurate  representations  of  the  continuum 


models  at  high  frequency. 


3  Beam  Model 

3.1  Euler  Bernoulli  Beam 


In  this  section  we  review  the  well-known  setup  for  the  cubic  finite  element  approximation  of 
the  Euler- Bernoulli  beam  equation.  Although  this  method  can  be  found  in  various  parts  in 
the  literature,  the  combination  of  finite  elements  in  space  with  time  derivatives  and  control 
forces  requires  piecing  it  together  from  a  couple  of  sources,  hence  it  is  induced  here  for  the 
sake  of  completeness. 

Consider  an  elastic  beam  of  length  l  with  uniform  cross-sectional  area  which  is  small 
compared  to  the  length.  The  distributed  force  acting  on  the  beam  is  assumed  to  be  zero. 
The  lateral  vibration  of  the  beam  can  be  described  by  Euler  Bernoulli  beam  equation: 


d*u>(t,x)  d2w(t,x)  ,dw(t,x) 

EI~d?~  +  p— S3-  +  k~dT'  ~  0  (1) 

where  0  <  x  <  /,  E,  I  axe  the  modulus  of  elasticity  and  the  moment  of  inertia  respectively, 
p  =  A  -  m  and  A  is  the  cross  sectional  area  of  beam  and  m  is  mass  per  unit  volume,  and  k 
is  the  viscous  damping  coefficient.  w(t,x)  is  the  vertical  deflection  of  beam. 

Several  cases  of  the  boundary  conditions  for  the  beam  are: 

a) ,  clamped  at  one  end  and  free  at  another,  ("clamped- free”),  controlled  by  a  shear  force 
at  the  free  end,  and/or  by  a  torque  at  the  free  end. 

b) .  hinged  at  both  ends,  ("hinged-hinged”),  controlled  by  a  torque  at  the  left  end  (x=0). 

c) .  hinged  at  one  end  and  free  at  another,  ("hinged- free”),  controlled  by  a  torque  at  left 
md  (x=0). 

When  there  is  no  external  force  acting  on  the  beam,  the  boundary  conditions  are  as 
follows: 


(i)  Free  end: 


Bl&jid  =  0,  «te  =  l 

ox1 
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(ii)  Hinged  end: 


(iii)  Clamped  end: 


d  .dw2(t,x) 

di{  dx 2 


)  =  0, 


atx  =  / 


w(t,  0)  =  0 


El 


dw2(t,x) 
dx 2 


=  0, 


atx  =  0 


tu(f,0)  =  0 


=  0,  atx  =  0 

3.2  Control  Inputs  and  Outputs 

We  assume  that  the  control  acts  on  the  beam  through  shear  force(s)  applied  at  a  point,  or  a 
bending  moment  (an  external  torque)  applied  to  one  or  both  ends  of  the  beam.  More  recently, 
distributed  force  actuators  made  out  of  piezo-ceramic  materials,  have  appeared.  Control  by 
using  such  actuators  can  be  modeled  by  an  easy  modification  of  the  finite  element  model 
described  below. 

Typical  configurations  of  control  actuators  and  sensors  are: 

A).  Clamped-free  beam  controlled  by  a  point  force  applied  to  the  free  end  (x  =  /). 
Boundary  conditions: 


dw(t,x) 

dx 


w(t,  0)  =  0 

§>.0  =  ° 


Control  input: 


Measured  outputs: 

1 

| 

yi(<)  =  w(ttl) 

«(<)  =  §;(<>') 

1 

B).  Hinged-hinged  beam  controlled  by  a  torque  applied  to  one  end 
Boundary  conditions: 

1 

■ 

to(i,0)  =  w(tj)  =  0 

-  • 

■ 

Control  input: 

1 

II? 

c*. 

*# 

o 

II 

<■» 

1 

Measured  outputs: 

1 

1 

1 

#.(<)  =  f|(<.0) 

«*>  - 

yz(t)  = 

»<«)  = 

1 

where  i,  is  the  sensor  location. 

C).  Hinged-free  beam  is  like  case  B)  except  the  condition  w(t,  /)  =  0 

is  replaced  by 

I 

d3w/dx3(t,l)  =  0. 

■ 

Let  the  state  space  for  the  system  be 

1 

(w(t,x),wt(t,x))  G  H2(0,l)  x  L3(0,l) 

(2) 

1 

■ 

II  II3- i<  a*1  )<ll  +  i(  «  ),ir 

(3) 

1 

Then  in  both  cases  A)  and  B)  the  control  input  operator  is  unbounded. 

1 

1 

7 

However,  in  case  the  control  torque  acts  on  the  beam’s  end  through  a  hub  with  a  nonzero 
moment  of  inertia  /«,  the  linearized  overall  system  has  a  bounded  input  operator: 


dw(ty0) 

dx 


—  4>(t)j  (*>  €  H*,wM  €  X2 


with  an  auxiliary  differential  equation 


Ililo~EI~d^~~u{t) 


(4) 


(5) 


where  is  the  angle  of  hub  rotation,  and  u(t)  is  a  control  torque. 


3.3  Kelvin- Voight  Damping  and  Rayleigh  Damping 

In  some  beam  models  the  damping  is  assumed  to  be  proportional  to  the  strain  rate,  thus 
contributing  the  term 


dx  2dt 


(6) 


where  c  is  the  Kelvin- Voight  damping  coefficient,  and  I  is  the  moment  of  inertia.  This 
term  is  easy  to  include  in  the  finite  element  model  described  below.  This  type  of  damping 
moves  the  high  frequency  eigenvalues  into  the  negative  real  semiaxis. 

The  linear  combination  of  (6)  and  viscous  damping  is  known  as  Rayleigh  damping. 


3.4  Finite  E  -ment  Model 

Finite  Element  Method  is  u?v*  to  approximate  the  continuum  Euler-Bernoulli  model  with 
a  finite-dimensional  model,  which  converts  the  spatial  and  time-  differential  equation  into 
time  differential  equation.  1 

The  beam  is  divided  into  N  segments  of  equal  length,  with  N  +  1  nodes,  as  shown 
in  Fig.(l)  At  any  arbitrarily  fixed  time  t  >  0,  the  displacement  over  each  element,  w,  is 
approximated  by  a  linear  combination  of  four  cubic  interpolation  functions  <f>j,  with  Uj  being 
the  value  of  w,  at  time  t,  at  both  ends  of  the  element. 


1  Detail  see  Appendix  A 


X 


^O3©  3  ...  •©•*'...  * 

Figure  1:  Finite  element  discretization  of  a  one-dimensional  beam. 


where 


+  rfM0  +  4e)4e) 


«1  =  10^  =  U>(Ze)  t*3  =  u4'*  =  tu(xe+l) 


«2  =  ^  =  &(Xe)  U<=0<e)  =  <?(*,+,)  (0,  =  -~ ) 


rie)  =  1-3(1-^)2+2(^t^)3 

he  /»e 

#  =  -(*-*,)(l-^)2 


4*’  =  3(£*7i),-2(f-*71)3 

4*’  =  -(*-:r.)[(^)J-^]  0) 

Using  the  approximation  of  the  form  of  Eq.(7)  for  w,  and  v  =  <f>j  in  the  variational  form 
of  Euler-Bernoulli  equation,  we  obtain  the  the  finite  element  model  over  element  e 


,x  -  x , 


M(e)u(e)  +  if  <e>u(e)  -1-  L<e>u(e)  =  Q(e) 
where  and  are  4  by  4  coefficient  matrices,  with 
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(11) 


Mjf  =  £'*'  pt^jix 

Kif  =  P'kMH* 


(e) 


-L 


r«  dx*  dx* 


and  is  the  boundary  condition  vector  for  element  e,  with 


Of* 

<??> 

<# 

q':' 


-fs<£0— 

-in¬ 


fill 


The  variables  tt^  in  Eq.(10)  represent  the  displacements  and  minus  slopes,  at  any  time 
t  >  0,  at  both  ends  of  element  e 


«w  =  [uie)t4e)t4e)u<e)]r 

=  [  w(xe)  0(xe)  w(ze+i)  0(ze+1)  ]T 


(13) 


Now  the  element  model  is  converted  into  a  second  order  time  differential  equation.  For 
the  case  in  which  El,  p  and  k  are  constant  over  an  element,  the  coefficient  matrices  M^e\K^ 
and  are  given  by  (Ke  corresponds  to  viscous  damping). 


'  156 

-22  h 

54 

13  h 

M(e)  = 

ph 

-22  h 

4  h2 

— 13h 

-3  h2 

420 

54 

-13  h 

156 

22  h 

.  13fc 

-3  h2 

22  h 

4  h2 

‘  156 

-22  h 

54 

13  h 

K{e)  = 

kh 

-22  h 

4  h2 

-13h 

-3  h2 

420 

54 

-13  h 

156 

22  h 

13h 

-3  h2 

22  h 

4  h2 

(14) 


(15) 
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■  6  -3 h  -6  -3 h 

(e)_2£/  -3  h  2  h2  Zk  k2 
L  ~  h*  -6  Zh  6  3 h 

.  -Zh  k 2  Zk  2 h2 

where  h  is  mesh  size. 

The  global  model  is  obtained  by  assembling  all  element  models  in  the  way  that  the  global 
coefficient  matrices  and  boundary  conditions  are  formed  by  overlapping  element  coefficient 
matrices  at  all  intermediate  nodes,  where  local  variables  are  repeated  for  two  consecutive 
elements,  such  as  tij?  =  =  uj,+1^  and  tt^  =  6^  =  =  Uj+1\  Hence,  the 

global  model  is  in  the  form 

MU  +  KU  +  LU  =  Q 

where  M ,  K  and  L  are  matrices  of  dimension  (2 N  -f  r),  (r 
(2 N  +  r)  vector.  The  global  variable  vector  U  is 

u  =  it/,  (/, . tw-iiw,  ,r 

=  [  <0(i|)  8(*i) . “’Otw+i)  ]T  (18) 


An  example  of  a  2  element  beam  is  shown  as  follow 


m<?  M#  Mj?  M{?  oo]  qS1’ 

m£>  M&>  M$?  0  0 

Mj?  M3(?  M&]  +  M,(?  Mff  +  Mj?  M<?  M<?  q_  qP  + 

m8>  M$  M$  +  m£]  m£>  +  M<?  M$  M<?  Q  Q{i]  +  Q?] 

0  o  m£>  q™ 

0  0  Ml?  Ml?  Ml?  Ml?  [  Q(42) 

Matrices  K  and  L  are  built  in  the  same  way  as  M.  We  obtain  the  global  model  for  the 
two-element  beam 

zIt  depends  on  boundary  condition 


(17) 

=  0, 1  or  2)a,  and  Q  is  a 
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Mh. 

420 


kh 

420 


‘  156 

-22  ft 

54 

13ft 

0 

0  ' 

ox 

— 22A 

4  ft2 

-13ft 

-3ft2 

0 

0 

u7 

54 

-13  ft 

312 

0 

54 

13ft 

o3 

13ft 

-3  ft2 

0 

8ft2 

-13ft 

-3ft2 

4 

u4 

0 

0 

54 

-13ft 

156 

22ft 

Us 

0 

0 

13ft 

-3ft2 

22ft 

4ft2 

.  Us  t 

'  156 

-22  ft 

54 

13ft 

0 

0  ‘ 

Ur 

—22ft 

4ft2 

-13ft 

-3ft2 

0 

0 

U3 

54 

—13ft 

312 

0 

54 

13ft 

u3 

13  ft 

-3ft2 

0 

8ft2 

-13ft 

-3ft2 

4 

U4 

0 

0 

54 

-13ft 

156 

22  ft 

Us 

0 

0 

13ft 

-3ft2 

22  ft 

4ft2 

.  Us  . 

2  El 
~W 


‘  6 

-3  ft 

-6 

-3ft 

0 

0  ' 

Ur  I 

'  F0  ' 

-3  ft 

2ft2 

3  ft 

ft2 

0 

0 

U3 

M0 

-6 

3ft 

12 

0 

-6 

-3ft 

u3 

0 

-3  ft 

ft2 

0 

4ft2 

3ft 

ft2 

< 

U4 

»  =  . 

0 

0 

0 

-6 

3ft 

6 

3  ft 

U5 

~Ft 

.  0 

0 

-3ft 

ft2 

3ft 

2ft2 

Us  J 

-Mi 

where 


F0  =  Of 
M„  =  of 
F,  =  Of 
M,  =  Of 


= 

-  -&«£»- 
=  -1^1- 


3.5  Imposition  of  Boundary  Conditions 

The  finite-element  model  in  Eq.(17)  is  valid  for  any  beam  described  by  the  differential  equa¬ 
tion  (in  Eq.(l)),  irrespective  of  the  boundary  conditions.  The  coefficient  matrix  (in  Eq.(17)) 
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is  singular  prior  to  the  imposition  of  the  essential  boundary  conditions  (on  the  primary 
variables).  Upon  the  imposition  of  suitable  boundary  conditions  of  beam,  we  obtain  a  non- 
singular  matrix.  First,  we  note  that  the  natural  boundary  conditions  are  included  in  the 
column  vector  {Q^}.  At  all  global  nodes  between  the  boundary  nodes,  the  sum  of  the  con¬ 
tributions  of  the  boundary  conditions  from  right  node  of  element  e  and  left  node  of  element 
e  -f  1  is  zero 


Q?  +  «S'+"  =  +  [EI  =  o  (20) 

The  boundary  conditions  at  two  ends  of  the  beam  are 


e!11  = 

Qi"  = 

Qf  = 

of*  - 

(21) 

For  hinged-free  beam,  the  only  known  essential  boundary  condition  is 


Ut  =  ti}1)  =  0  (22) 

Using  Eqs.(20)  and  (22)  in  Eq.(17),  we  have 


’  {A-11} 

{a:12}  ‘ 

'  &  =  o ' 

{ K 21) 

{ K 22} 

• 

• 

• 

. 

.  U2N+7 . 

13 


r  {l") 


{£,»} 


{L»}  1 


Ui  —  O 


{Ln) 


U2 

U2N+2 


F0 


Mo 

0 

• 

0 


-Ft  =  0 
-M,  =  0 


(23) 


where  each  coefficient  matrix  is  divided  into  four  parts  corresponding  to  variable  partition. 
The  fact  that  the  first  variable  U\  is  known  and  the  remaining  (2 N  +  1)  variables  are  to  be 
determined  provides  us  the  motivation  to  partition  (shown  by  dashed  lines)  the  matrix 
equation  (23).  Eq.(23)  is  in  the  form  (where  V1  =  U\). 


'  {A/u} 

)M'2) 

/  V*  \ 

IK'2)  1 

1  Y1  \ 

.  f {£“} 

{M21} 

(  M22} 

\  V2  ) 

[  {^21} 

{K22} 

l  V2/ 

4  {£’■} 

Equation  (24)  can  be  written,  after  carrying  out  the  matrix  multiplication,  in  the  form 


{mu}wi  +  {Af12}v2  +  {Arn}vi  +  {K12}v2  +  +  {z,12>v2  =  { q *}  (25) 

{M21}V*  +  {M22}V2  +  {K2'}V 1  +  (tf22}V2  +  {L21}?1  +  {£22}V2  =  {Q2}  (26) 

Since  {V1}  and  {£2}  are  known,  we  can  use  Eq.(26)  to  solve  for  {V2}.  In  other  word, 
we  can  simply  cross  out  the  first  row  and  the  first  column  of  the  coefficient  matrices  in  this 
case  to  get  the  modified  model. 

For  hinged-hinged  beam,  the  deflections  at  both  ends  are  zero 

Ux  =  u<11)  =  0 

U2N+ 1  =  uf  >  =  0  (27) 

Using  same  procedure,  we  obtain  modified  model  for  hinged-hinged  beam  by  deleting  the 
1st  and  (2 N  +  l)th  rows  and  columns  from  coefficient  matrices. 
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For  clamped-free  beam,  the  deflection  and  slop  at  left  end  are  zero 

Ut  =  u(,x)  =  0 

U2  =  u[a)  =  0 

(28) 

The  modified  model  is  obtained  by  deleting  the  first  two  rows  and  columns  from  coefficient 
matrices. 

For  clamped-clamped  beam,  the  deflections  and  slopes  at  both  ends  are  zero 


Ui 

II 

c 

** 

II 

© 

u2 

© 

II 

S- 

3 

II 

UtN+ 1 

=  4N)  =  o 

UiN+i 

o 

II 

3 

II 

(29) 

The  modified  model  is  obtained  by  deleting  the  first  two  rows  and  columns,  and  the  last 
two  rows  and  columns  from  coefficient  matrices. 


3.6  Finite  Dimensional  Control  Model 

After  incorporating  the  boundary  conditions  and  output  equations,  we  obtain  the  following 
finite  dimensional  control  model 


Mz  +  Kz  +  Lz  =  Gu(t) 

J,(0  =  iST  f  ^  1 


(30) 

(31) 


where  z(t)  =  V3(t),u(t)  =  control,  y(t)  —  output.  (Note  that  u(t)  is  different  from  vector 
U  in  equation  17).  Matrices  M,K ,  and  L  are  reduced  versions  of  matrices  M,  K,  L  appearing 
in  (17),  the  reduction  depending  on  the  boundary  conditions. 

Case  A).  For  a  clamped-free  beam  with  N  finite  elements,  the  size  of  matrices  Af,  Ky  L 
is  2 N  x  2  N. 
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G  is  a  2N  x  1  column  vector  which  is  zeros  except  for  G(2N  —  1)  =  1. 

If  y(t)  =  w(/,  t),  then  H  is  a  1  x  4 N  row  vector  which  is  zero  except  for  H(iN  —  1)  =  1. 
Case  B).  For  a  hinged-hinged  beam  with  control  torque  at  the  left  end  (x=0)  of  the  beam, 
and  measured  time  derivative  of  the  beam’s  slope  at  x=0,  the  size  of  matrices  M,  K ,  L  is 
2  N  x  2  N. 

G  is  a  2 N  x  1  column  vector  with  (7(1)  =  1  and  G(j)  =  0  for  j  >  1. 

H  is  a  1  x  4 N  row  vector  with  H(2N  +  1)  =  1  and  all  other  entries  equal  to  zero. 
Controllability  of  the  finite  dimensional  models  can  be  verified  by  a  modified  Hautus 
condition 


Rank  [MX2  +  KX  +  L\G]  =  2 N,  VA  6  C  (32) 

For  a  beam  with  zero  damping,  the  condition  simplifies 

Rank  [-Mu2  +  L\G\  =  2 N,  VuZK  (33) 

Numerical  tests  performed  in  the  examples  investigated  in  this  report  showed  this  con¬ 
dition  to  be  always  satisfied  in  Case  A  and  Case  B. 

4  The  DIRK  Method  for  the  Integration  of  Finite 
Element  Model 

In  the  process  of  numerical  integration  of  the  beam  equations,  one  needs  to  use  schemes  that 
are  suitable  for  highly  oscillatory  systems.  In  addition,  one  needs  to  accommodate  the  fact 
that  the  control  input  may  be  in  the  feedback  form  rather  than  an  exogenous  time  function. 
That  is,  at  a  given  time  t  only  the  current  value  of  control  u(t)  is  available,  while  the  future 
values  ( u(t  -f  <5),<5  >  0)  of  control  have  not  yet  been  determined  by  the  control  algorithm. 

The  finite  element  model  described  in  section  3.4  is  a  system  of  second  order  differential 
equations.  This  system  can  be  transformed  into  a  standard  state-space  form 

x  =  Ax  +  Bu  (34) 
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where 


A  ~  -Af~xL  0  B~  M~XQ  ^ 

In  numerical  calculation,  direct  use  of  matrices  M~XL  and  M~XQ  is  impractical.  It  is 
better  to  perform  a  preliminary  transformation  of  variables  that  maintains  the  symmetry  of 
submatrices  appearing  in  matrix  A.  Let  M  =  McMj  denote  the  Cholesky  decomposition  of 
matrix  M.  Then  in  reference  to  equation  (17),  let 


MjV 


m;'Q 


Then  the  original  system  is  replaced  by 


y  +  Kmy  +  Lmy  =  Qm 

In  this  case  the  corresponding  first  order  system  in  state  space  form  has  matrices 


The  matrix  A  is,  for  the  Euler  Bernoulli  equation,  a  poorly  conditioned  matrix.  For  ex¬ 
ample,  a  finite  element  model  with  4  elements  will  have  matrix  A  with  cond  (A)  =  1.2902e+4. 
For  64  elements  the  corresponding  number  is  cond  (A)  =  8.4557e+8.  As  a  result  of  this,  the 
numerical  accuracy  of  some  operations  involving  A  decreases  with  the  increase  in  the  num¬ 
ber  of  elements,  making  the  numerical  computation  difficult  for  high  number  of  elements.  A 
small  step  size  is  required  to  keep  track  of  high  frequencies,  but  a  small  step  size  increases 
the  round-off  errors.  As  a  result  of  this,  the  use  of  standard  numerical  integration  methods 
for  the  state  equations  leads  to  very  large  errors  (e.g.  the  Matlab  subroutine  ’lsim’  fails  for 


N  >  10).  However,  a  recently  developed  extension  of  the  Diagonally  Implicit  Runge  Kutta 
(DIRK)  method  to  second-order  systems  exhibits  a  very  good  numerical  accuracy.  The  key 
point  here  is  that  this  forces  the  use  of  the  second  order  systems  rather  that  the  augmented 
first  order  system. 

4.1  DIRK  Method  for  Oscillatory  Second-order  Differential  Equa¬ 
tions 

The  DIRK  method  proposed  by  Van  der  Houwen  and  Sommeijer  (26)  is  devised  for  second- 
order  differential  equation  without  damping 

y  =  f(t,y ),  y(0)  =  yo,  y(0)  =  yo  (39) 

The  general  m-stage  DIRK  method  for  the  system  of  ODEs  (Eq.(39))  is  given  by 

m 

Yni  =  yn  +  CjAyn  +  h^aji/^n  +  ciA,!^,),  j  =  l,...,m 

/= 1 

yn+i  =  yn  +  hyn  +  h2  51  bjf{tn  -f  Cjh,  Ynj) 

i=i 

ifn+l  =  yn  +  h^2  *>'/(<„  +  Cjh,  Ynj)  (40) 

3= 1 

where  the  parameters  Ojj,  6j,  and  Cj  are  assumed  be  real.  By  defining 

A:=(«ii),  6:=(6i),  ,6':=  (6'),  c  :=  (Cj)  (41) 

the  DIRK  method  can  be  represented  by  the  Butcher  array: 

c  A 

P"  (42) 

V1 

An  example  of  Butcher  array  for  two-stage  DIRK  method  with  algebraic  order  p  =  2  is 
shown  as 
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The  system  with  damping  is  described  by 

y  =  /(*»y,y)i  y(0)  =  yo,  y(0)  =  y0  (44) 

A  DIRK  method  for  such  system  is  not  found  in  present  literature.  By  extending  the 
results  of  Van  der  Houwen  and  Sommeijer  for  zero-damping  system,  we  obtain  the  DIRK 
method  for  system  with  damping 

m 

Ynj  =  yn  +  Cjhyn  +  £  a]tf(tn  +  c,h,  Yni ,  Ynl),  j  =  1, m 

{=1 

m 

Ynj  ~  jin  d~  ^  ^  ,  Qjlfjtn  Cjh,  Ynl,  Ynl ),  J  = 

1=1 

m 

yn+1  =  y«  +  Ayn  +  A2  £  4j/(/n  +  Cjh,  Ynj ,  v„/) 

J=1 

yn+1  =  J/n  +  hf:6'/(tn  +  cAKi,Knl)  (45) 

;=i 

where  the  parameters  aji ,  bj,  6'  and  Cj  are  from  Butcher  array,  and  a',  are  chosen  as  same 
as  aji.  The  simulation  results  show  that  the  extended  DIRK  method  works  with  reasonably 
good  accuracy. 

4.2  DIRK  Formulas  for  the  Beam  Problem 

The  finite-element  model  (Eq.(17))  for  Euler  Bernoulli  beam  can  be  written  as  following 

y  =  /(y,y) 

=  -M-'Ky-M-'Ly  +  M-'Q 

-  ~Kmy  -  Lmy  +  Qm  (46) 
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where  Km  —  — M~lK ,  Lm  =  M~XL  and  Qm  =  M~XQ.  Actually,  to  preserve  the  symme¬ 
try  of  the  matrices,  we  later  use  Cholesky’s  decomposition  before  the  inversion  of  M. 

Using  DIRK  method  (Eq.(40))  with  m  =  2,  we  obtain  the  discrete  time  model  for  the 
system  without  damping  (K  [0]) 


Yni  =  y*»  +  c\hyn  4-  /i2[aji(— LmYni  +  Qm)  +  ai2(— LmYni  4-  C?m)]  (47) 

Yn 2  =  yn  +  c2hyn  4-  A2[a2i(— LmYn\  +  Qm)  4-  a^— LmYn2  -f  Qm)]  (48) 

J/n+l  =  Vn  +  hyn  +  A2[6i(  —  LmYn\  +  Qm)  +  6s(  —  LmYn2  4-  Qm)\  (49) 

y„+i  =  y„  4-  h[b\(-LmYnl  +  Qm)  +  b'2(—LmYn2  4-  <?m)]  (50) 


The  discrete  time  model  (Eq.(49)  and  Eq.(50))  still  have  variables  Yni  and  Yn 2  in  their 
expressions.  To  simplify  the  model,  we  solve  Eq.(47)  and  Eq.(48)  to  get  Ki  and  Yn2.  Noting 
the  Yn  1  and  Yn 2  are  also  functions  of  yn  and  yn  and  au  =  0,  we  express  the  Yn  1  and  Yn2  in 
the  form 

Yn  1  =  A\yn  4-  A2yn  4-  AnQm 

Yn2  =  A$yn  4-  A4yn  -f  A22Qm  (51) 

where 

A,  = 

A2  =  Cih(I  +  h'2a11Lm)  1 
An  —  h2an(I  +  h2anLm)  1 
A3  =  (/  +  h2a22Lm)~l(I  -  h^LnAj 
A+  =  (/  4-  h2a22im)_1(c2h/  —  h2a2\LmA2) 

Aj2  —  (I  +  h2a22Lm)~1(—h2a2iLmAu  4-  h2{a2\  -4  <*22)/)  (52) 
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and  /  is  identity  matrix  with  same  size  as  Lm. 

Substituting  Yn i  and  Yn 2  in  Eq.(49)  and  Eq.(50),  we  obtain  the  following  recursive, 
discrete-time  model 

yn+i  =  yn  +  hyn  -  h2[biLm(Aiyn  +  A2yn  +  AnQm) 

+  b2Lm{.Azyn  +  A4yn  +  A22Qm)  +  (&i  +  6a)Qm] 

=  [/  -  h’MMi  +  M3)]y»  +  [hi  -  h2Lm(bxA2  +  M<)]y„ 

+  [h2(6i  +  bi)I  —  h2  Lm(biAu  +  iMaaJJQm 

J/n+1  =  yn  -  h[b\Lm(Aiyn  +  A2yn  +  AnQm) 

+  ^Lm^A^yn  +  A4yn  +  A22Qm)  +  (^1  +  ^2)Qm] 

=  -hLm(b[Ai  +  b’2A3)yn  +  [/  -  hLm(b\A2  +  6'A<)]yn 
+  [h(b\  +  h'7)I  -  hLm(b\Au  +  6'2A22)]Qm  (53) 

Using  the  standard  state-space  form,  we  have 


z„+i  =  Fx  „  +  Gun  (54) 

zn  =  Hxn  (55) 


where  xn  =  [y„  yn]T  is  a  (4 N  +  2r)  by  1  vector,  standing  for  displacements  and  minus 
slopes  of  all  element  (see  Eq.  (17))  and  their  velocities,  u  is  input  as  bending  moment,  and 
z  is  output  which  is  measurable, 
and 


/  -  fc’MMx  +  Ms)  hi  -  A2Lm(M3  +  M<)  ‘ 
-hLm{b\Ax  +  b'2A3)  I-hLnMAt  +  b^At) 


P  __  h2(bi  +  hi)I  —  h2 Lm{b\A\\  +  6a/l22]^m 
~  h{Vx  +  V2)I-hLm{b\Au  +  V2A22)\Qm  . 

(57) 

For  damped  system,  the  extended  DIRK  method  (Eq.(45))  is  used  to  get  following  equa¬ 
tions 

Yn  i  =  yn  +  c\hyn  +  h2[an(—KmYni  —  LmYni  +  Qm)  +  <*i  i{—KmYn2  —  LmYn2  +  Qm)] 
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Ynt  =  3/n  +  c2/iy„  +  /x2[a2i(— /^m^ni  —  LmYni  +  Qm)  -f  a2j(—KmYn2  —  LmYnj  +  Qw)] 
^nl  =  J/n  +  MaXl(~  KmYnl  ~  LmYnX  -f  Qm)  +  <*xj(~^i»*Vn2  ~  LmYn2  +  Qm)] 

K»2  —  J/n  +  Maji(-"-K”»Vni  ~  LmYni  -4*  Qm)  +  o27(—KmYn2  ~  LmYni  +  Qm)]  (58) 


I/n+1  —  Vn  +  h.yn  +  h  [6i(  KmYn\  —  LmYn\  +  Qm)  ■+■  bl^—KmYn'l  —  LmYn2  +  Qm)] 
iln+l  =  j)n  +  M^l(  —  KmYn\  —  LmYn  1  +  Qm)  +  b2(—KmYn2  —  LmYn2  +  Qm)]  (59) 

Solving  equations  in  Eq.(58),  we  have  KnJ,  Yn 2,  K„i  and  YnX  in  the  form 


where 
A-ki  = 

£*x  = 


[£’]=‘4“[£]+Bm<?" 


I  +  h2ai\Lm 

h2axlKm 

-1 

J  cxhl 

^all  Lm 

I  +  ka’nKm  . 

0  I 

I  +  h^aiiLm 

h2auKm 

-1 

A2Oll I  ] 

ha’nLm 

I  +  ha'nKm 

ha'^I  J 

(60) 

(61) 


Al(2  — 


Bk  2  = 


/  +  ha’llKm  h2a.j2Km 
^  “h  ^®2J^m 

/  +  ha'^Krn  h2a,22Km 

ha^Lm  I  “I”  ^®22^m 


1-1  r 


-1 


/  cl/lJ 
0 


A/  —  >* 


h7a2\Lm  h2a2iK„ 


ha'uK„ 


Am) 


h2(a  2i  +  022)/ 
h(a21  +  022)-^ 


h2a2iLm  h?a2iKm 
ha'21Lm  ha2lKm 


Bkl) 


(62) 


Substituting  Eq.(60)  and  Eq.(61)  into  Eq.(59)  and  rearranging  leads 
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/  W  1  r  ».  I  _  f  h'b,L »  h%K„  1  f  Ki  1 

■  l  0  /  J  l  in  J  l  k%Lm  hHKm  J  l  n,  J 

'  k%Lm  h%Km  I  f  Y.7  1  f  M(4,  +  6,)/  1  - 

hV2U  h%K„  J  l  Yn,  J  +  l  h(Vt  +  V,)I  J  V" 


_  r  /  w  i  r ».  i  r  ***,*«  i , .  s„  i ,  „  0  , 

-  0  /  J  l  in  J  -  [  kHL.  kUK.  j  (A*‘  [  i,  J  +B“<?”) 


'  h2b 
hb'2 


h2kiLm  h*biK„ 
hb'2Lm  hb2Km 


Vn 

.  y » . 


+  BkiQm)  + 


h\b X 


bi  +  6j)/  ' 

{  +  *4)/  . 


[  /  hi  I  f  /i2&iAr, 

~  0  /  /liiln,  A6- 


Wrn]A  _  h2b2Lm  h2b*Km  1  [  yn  ' 

\Km  Akl  hb'2Lm  hbf2Km  Ak2)  yn 


Aa(*i  +  *a)/  1  f  h% Lm  h%Km  ]  _  r  h*biLm  h'biK, 


h(K  +  b'2)I 


hb\Km  I  I  M&Im  hV2K, 


Wm  ' 
\Km 


Bki)Qr 


In  the  standard  state-space  form,  we  have 


where 


I  hi 
0  I 


_/[  Wlh 
~([  h(b i 


+ 

+  *»/ 


*n+l  =  FXn  +  Gun 

Zn  ~  H  %n 


h%Lm  h'hKrn  In  _  [  h*b2Lm  h'biKn 
hb\Lm  hV2Km  Akl  hb2Lm  hV2Km  Ak2 


hKm 

V2Km 


h%Lm  h%K , 
hb[Lm  hb\K„ 


Km  d  hPbiLm  bPbiKm  „  ((.R. 

\Km  \  Bk'  ~  [  hV2Lm  WtKm  J  Bu)Q -  (66) 


5  Numerical  Results  of  Simulation 


Numerical  simulations  have  been  done  for  several  cases,  such  as  initially  excited  mode,  feed¬ 
back  control  and  observer  design.  In  each  case,  various  beams  with  different  boundary  con¬ 
ditions  were  used,  which  are  hinged-hinged,  hinged-free,  clamped-free  and  clamped-clamped 
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beam.  The  finite-element  model  and  the  DIRK  method,  described  in  section  3  and  4  ,  were 
used  in  the  simulations.  Also,  some  comparisons  between  numerical  results  and  analytical 
solution  have  been  done. 

5.1  Beam  Parameters 

The  beam  used  in  simulation  is  Euler-Bernoulli  beam  (Eq.(l))  with  a  finite-element  model. 
Physically,  the  following  thin  aluminum  alloy  beam  corresponds  to  the  combination  of  pa¬ 
rameters  used  in  our  simulations 
Length  /  =  lm 

Cross-section  width  wc  =  6  x  10_3m 
Cross-section  height  hc  =  3.01  x  10~3m 
(the  beam  vibrates  in  the  "vertical”  direction) 

Cross  sectional  area  A  —  1.809  x  10“6m2 
Area’s  moment  of  inertia  /  =  1.37  x  10~nm4 
Young’s  modulus  E  =  7.309  x  10 7kg/m.sect 
Mass  density  m  =  2.768  x  103A:<//m3 
Then  p  =  mA  =  5.0072  x  10 2kg/m 
El  —  1.0015  x  10 ~3kg.m3/sec2 
El  Ip  =  0.02  m4/sec2 

The  0.02  ratio  corresponds  to  a  thin,  fairly  elastic  beam,  which  exhibits  several  vibration 
modes.  We  could  have,  of  course,  chosen  a  different  value  of  El/ p,  which  would  yield  a 
less  elastic,  and  thus  easier-to-control  beam,  but  then  some  of  the  essential  features  of  such 
control  problem  would  be  lost. 

5.2  Analytical  Solution  of  Euler  Bernoulli  Beam 

The  undamped  Euler  Bernoulli  beam  equation  is 


2^u;(<,x)  03u>(*,x) 

c  ~a?~  +  -*r-  - 0 
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where 


c  = 


(68) 


The  analytical  solution  of  Eq.(67)  with  constant  boundary  conditions  can  be  expressed 
in  Fourier  series  form: 


^(r,t)=E^«(*)rn(<)  (69) 

n=t 

The  functions  Tn(t)  can  be  expressed  as 

Tn(t)  =  Ancosunt  +  Bnsinunt  (70) 

where  An  and  Bn  are  constants  that  can  be  found  from  the  initial  conditions. 

The  modal  shape  function  Wn(t)  are 

Wn(x)  =  Cicosfinx  +  Cisin/3nx  +  C3coshf)nx  +  CAsinhpnx  (71) 

where  Ci,  C3,  C3,  and  C4  are  constants  that  depend  on  the  boundary  conditions.  The 
natural  frequencies  of  the  beam  are 


The  modal  shape  functions  and  values  of  0„l ,  corresponding  to  different  boundary  con¬ 
ditions,  are  as  follows 

(i)  Hinged-hinged  Beam 


Wn(x)  =  sinfinx 


(73) 


PJ  =  (ri  +  l)v,  n  =  0,1,2,... 


(ii)  Clamped-free  Beam 


(74) 
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W,(x)  =  coshfa  -  CM/*.*  -  + 


(75) 


/*,/  =  n*  +  5 +«<-).  "  =  0,1,2,...  (76) 

Z  n 

5.3  Comparison  of  natural  and  finite-element  modal  frequencies 

In  all  the  numerical  comparisons  given  below,  the  numerical  parameters  of  the  beam  were 
/  =  1  m,EI/p  =  0.02,  Ar(damping)=  0. Interpretation  of  this  choice  is  given  in  section  5.1. 

The  tables  below  show  eigenfrequencies  u>4,wg,fa;i6,W32ju>64  obtained  from  finite  element 
models  with  4,  8,  16,  32  and  64  elements  respectively,  compared  with  u>n  =  frequencies 
corresponding  to  the  PDE  model. 

These  tables  confirm  the  fact  that  the  first  few  natural  frequencies  are  modeled  by  FEs 
quite  accurately,  while  the  high  frequencies  are  not.  The  number  of  accurately  modeled 
frequencies  increases  with  the  number  of  finite  elements. 


Mode  No 

U)A 

W8 

Wl6 

^32 

^64 

“>»  1 

1 

1.3961 

1.3958 

■ 

1.3958 

1.3958 

2 

5.6051 

5.5845 

5.5831 

5.5831 

3 

12.7915 

12.5781 

12.5630 

12.5620 

12.5620 

12.5623 

4 

24.7871 

22.4205 

22.3382 

22.3327 

22.3324 

22.3324 

5 

39.3990 

35.2178 

34.9162 

34.8957 

34.8944 

34.8943 

6 

62.3044 

51.1660 

50.3125 

50.2520 

50.2481 

50.2478 

7 

93.3411 

70.5230 

68.5536 

68.4033 

68.3935 

68.3929 

8 

113.5887 

99.1484 

89.6820 

89.3527 

89.3309 

89.3295 

9 

- 

123.4418 

113.7595 

113.1044 

113.0606 

113.0576 

10 

- 

157.5959 

140.8713 

139.6649 

139.5829 

139.5773 

11 

- 

199.0853 

171.1279 

169.0429 

168.8984 

168.8885 

12 

- 

249.2176 

204.6639 

201.2499 

201.0079 

200.9913 

13 

- 

308.4743 

241.6236 

236.3008 

235.9125 

235.8856 

14 

- 

373.3646 

282.0920 

274.2146 

273.6133 

273.5715 

15 

- 

430.3208 

325.6174 

315.0143 

314.1120 

314.0489 

16 

- 

454.3549 

396.5935 

358.7281 

357.4106 

357.3178 

Table  1:  Comparison  of  natural  frequency  and  FE  modal  frequency  for  hinged-hinged  beam 


Mode  No 

re< 

re s 

»*ej6 

re3 7 

rc«4  i 

1 

2.59 66e-M 

1.6442e-ttS 

1. 0310c-1* 

“1.5899c-04  J 

2 

3.9469c-03 

2.5966e-04 

1.6443e-05 

1.0311c-06 

6.4215c-04  I 

3 

1.8273c"02 

1.2866c-03 

8.2786c"05 

5.2125c-06 

3.2625c"04 

4 

1.0992c"01 

3.9469c-03 

2.5966c-04 

1.6443c-05 

1.0311c-04 

5 

1.2909e-01 

9.2709c-03 

6.2781c-04 

4.0044c-05 

2.5156c"04 

6 

2.3994e-01 

1.8273c'02 

1.2866c-03 

8.2786c-05 

5.2125c"04 

7 

3.6478e-01 

3.1145e-02 

2.3508c"03 

1.5283c'04 

9.6482c-04 

8 

2.7157e-01 

1.0992e~O1 

3.9469c"03 

2.5966c-04 

1.6443c-04 

9 

- 

9.1849c"02 

6.2087c-03 

4.1402c'04 

2.6307c*04 

10 

- 

1.2909e_O1 

9.2709c-03 

6.2781c-04 

4.0044c-04 

11 

- 

1.7880e~O1 

1.3260c*02 

9.1401c-04 

5.8544c-04 

12 

- 

2.3994c"01 

1.8273c-02 

1.2866c-03 

8.2786c"04 

13 

- 

3.0773e"01 

2.4325e-02 

1.7603c-03 

1.1383c*04 

14 

- 

3.6478c-01 

3.1145e-02 

2.3508c'03 

1.5283c-04 

15 

- 

3.7023e"01 

3.6837e-02 

3.0741c-03 

2.0100c-04 

16 

- 

2.7157c"01 

1.0993c"01 

3.9469c"03 

2.5966c-04 

Table  2:  Relative  errors  of  frequencies  of  FE  models  for  hinged-hinged  beam 


Mode  No 

L)& 

<*>16 

<*>32 

<*>64 

<*>» 

1 

0.4973 

0.4972 

0.4972 

0.4972 

0.4972 

0.4972 

2 

3.1198 

3.1164 

3.1162 

3.1161 

3.1161 

3.1161 

3 

8.7929 

8.7306 

8.7257 

8.7253 

8.7253 

8.7253 

4 

17.3464 

17.1364 

17.1007 

17.0983 

17.0981 

17.0981 

!  5 

32.2635 

28.4279 

28.2759 

28.2651 

28.2645 

28.2644 

6 

51.8153 

42.7241 

42.2600 

42.2246 

42.2223 

42.2221 

7 

82.1445 

60.1431 

59.0730 

58.9781 

58.9718 

58.9714 

8 

134.7818 

79.8992 

78.7473 

78.5279 

78.5132 

78.5122 

9 

- 

112.3348 

101.3318 

100.8777 

100.8467 

100.8446 

10 

- 

142.4877 

126.8937 

126.0326 

125.9726 

125.9685 

11 

- 

180.7453 

155.5193 

154.0001 

153.8914 

153.8840 

12 

- 

227.8248 

187.3077 

184.7900 

184.6038 

184.5910 

13 

- 

285.2434 

222.3471 

218.4154 

218.1107 

218.0895 

14 

- 

352.5099 

260.6273 

254.8927 

254.4132 

254.3796 

15 

- 

420.0845 

301.6366 

294.2423 

293.5127 

293.4612 

16 

- 

541.9361 

341.5518 

336.4892 

335.4110 

335.3344 

Table  3:  Comparison  of  natural  frequency  and  FE  modal  frequency  for  clamped-free  beam 


Mode  No 


re4 


3.278c 

1.165c 

7.742c 

1.452c 

1.415c 

2.272e" 

3.930c 

7.167c 


re8 

rcie 

re  37 

re«4 

2.156c-06 

2.042c'07 

FWttJ 

1.372e'** 

8.002c-05 

5.165c'06 

3.776c -07 

7.439c"06 

6.084e-M 

3.987c'05 

2.623c'06 

2.689c"07 

2.240c'03 

1.511c'04 

9.574c"06 

5.559c"07 

5.786c-03 

4.084c'04 

2.642c'05 

1.855c"06 

1.189e-M 

8.979c'04 

5.829c'05 

3.671c"06 

1.987c'02 

1.723c'03 

1.133c"04 

7.168c"06 

1.767c'02 

2.995c'03 

1.998c"04 

1.269e"05 

1.140e-°l 

4.831c'03 

3.280c"04 

2.091c"05 

1.311c'01 

7.345c'03 

5.089c"04 

3.258c"05 

1.746e'01 

1.063e-02 

7.547c'04 

4.854e"05 

2.342e-01 

1.472e'02 

1.079c'03 

6.973c"05 

3.079e'°l 

1.952e-02 

1.494c"03 

9.716c"05 

3.858e'01 

2.456e-02 

2.017c'03 

1.319c"04 

4.315e'01 

2.786e-02 

2.662e-03 

1.752c"04 

6.161c'01 

1.854e-02 

3.444c-03 

2.283c"04 

Table  4:  Relative  errors  of  frequencies  of  FE  models  for  clamped-free  beam. 


5.4  Initially  Excited  Modes 

In  this  simulation,  the  behaviors  of  beam  are  examined  when  their  first  five  modes  are  excited 
separately  by  initial  functions.  Also,  the  accuracy  of  finite-element  model  and  DIRK  method 
used  in  Euler  Bernoulli  beam  problem  is  studied. 

Undamped  beams  with  two  types  of  boundary  conditions,  hinged-hinged  and  clamped- 
free,  are  considered  here.  A  beam  is  initially  bent  in  the  shape  described  by  initial  function, 
and  released.  The  lateral  vibration  of  beam  is  simulated  numerically,  and  the  reconstruction 
of  the  shape  of  beam  is  done  by  using  the  finite-element  model. 

5.4.1  Hinged-hinged  Beam 

In  the  simulation,  the  modal  shape  functions  were  used  as  initial  functions 

Wn(x)  =  Cn[$in/}nx\  (77) 

where  finl  =  (n  +  1)jt,  n  =  0, 1, 2, ...  .  With  C„  =  1,  we  have  initial  functions  for  the  first 
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Figure  3:  Vertical  motion  of  hinged-hinged  beam,  mode  3 
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the  first  mode.  Let  T\  and  T* 443  be  the  positions  on  the  {-axis  where  the  exact  solution 
assumes  its  first  and  4443th  zeros,  respectively.  The  time  interval  for  this,  T  =  Tu* 3  —  Tj, 
is  the  length  of  2221  oscillations  of  the  true  solution.  This  distance  is  compared  with  its 
numerical  analogue.  In  numerical  simulation,  the  distance  between  the  first  and  4443th 
"numerical  zeros”  is  denoted  by  T  —  T4443  —  T\.  The  period  preservation  of  DIRK  method 
in  the  beam  problem  is  examined  by  checking  how  close  the  T  is  to  T.  The  simulation  has 
been  done  for  4,  8,  16, and  32  element  beam  models  for  the  first  mode.  The  "numerical 
zero”  is  obtained  by  using  linear  interpolation  based  on  neighboring  numerical  values.  The 
simulation  results  are  shown  in  Table  5,  which  is  for  the  vertical  vibration  of  displacement 
in  the  middle  point  of  beam. 


In 

T  (sec) 

T(sec) 

T  -  f(sec) 

ID 

9. 9954 1 725840972 1  e+“ 

9  995422709364795e+O3 

-5 . 45095507732 1 851 e-*3 

8 

9.9978483013072096+°* 

9.997851666551 194e+03 

-3.365243985172128c-03 

16 

9.998002382497589e+03 

9.997997585413077e+O3 

4.797084513484151e-b* 

32 

9.998012037613589e+fl3 

9.998007922038640e+O3 

4.115574949537404c-03 

Table  5:  Period  Preservation  of  DIRK  method  for  hinged-hinged  beam 


The  energy  preservation  is  examined  by  comparing  the  norm  of  initial  state  variables 
(displacements  and  slopes  at  all  nodes,  and  their  velocities)  and  the  norm  of  state  variables 
after  2221  oscillations.  If  the  ratio  of  the  two  norms  equal  to  one,  the  DIRK  method  is  said 
to  preserve  energy  in  the  beam  problem.  Table  6  shows  the  norms  of  state  variables  at  initial 
time  and  the  2221th  periods.  The  simulations  were  done  based  on  4,  8,  16,  and  32  element 
beam  models  for  the  first  mode. 


nr 

Norm  1 

Norm  2 

Ratio  (n2/nj) 

■dr*>JVJUr*  mviMSium 

5.622207813156914e+°° 

1.0000064458247556+“ 

8 

7.3039730288006 1 2e+°° 

7^30399 177441 3425e+°° 

16 

9.840042663007326e+“ 

9.840095651 102425e+O(J 

1.000005384945667c+°° 

32 

1.355666901633727e+61 

l~355672308055335e+O1 

1. 00000398801 6232e+“ 

Table  6:  Energy  Preservation  of  DIRK  method  for  hinged-hinged  beam 
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Prom  the  simulation  results  (Table  5  and  Table  6),  we  can  see  that  the  DIRK  method 
used  in  the  finite-element  model  of  hinged-hinged  beam  preserves  both  period  and  energy  of 
oscillation  quite  well,  and  does  not  have  big  difference  for  models  with  various  sizes.  As  an 
example,  the  first  and  last  periods  of  trajectory  of  beam,  displacement  at  middle  point,  are 
shown  in  Fig.(5)  and  (6),  which  was  done  for  a  32  element  FE  model. 

The  accuracy  of  finite-element  model  plus  DIRK  method  is  measured  by  calculating  the 
£2  norm  and  the  energy  norm  of  the  error  between  exact  solution  and  numerical  solution 
which  are  based  on  FE  model.  The  two  norms  are  defined  by 

£2  norm: 

||e||o  =  Iitu«  -  t£)nu||o  as  {jf  [wtx  -  u>„„j2dx}1/2  (79) 

"Energy”  norm3 

IMU  =  IK.  -  «wlU  =  (/  -  %l  V*)1'1  (so) 

where  wex  and  wnu  are  exact  solution  and  numerical  solution  respectively. 

In  our  experiment,  the  "Energy  norm”  was  calculated  corresponding  to  each  mode. 
Hence,  the  exact  solution  u;ex  in  Eq.(80)  is  replaced  by  modal  solution  un 

wn  =  sinPnzcosu>nt  (81) 

The  numerical  solution  wnu  is  obtained  by  reconstructing  beam  shape  through  finite 
element  model. 

The  £2  errors  of  finite-element  model  and  DIRK  integration  method  for  the  first  five 
modes  are  shown  in  Table  7,  where  t\  to  U  are  1,  5,  8,  and  10  seconds  respectively.  The 
simulations  were  done  based  on  4,  8,  16,  32,  and  64  element  models,  using  different  step 
sizes,  baring  with  modes.  The  logarithm  of  the  £2  error  versus  the  logarithm  of  element 
number  at  t  =  Is  is  shown  in  Fig.(7),  and  the  logarithm  of  the  £2  error  versus  the  logarithm 
of  time  are  shown  in  Fig. (8).  The  simulation  results  show  that  the  finite  element  model  and 

3This  energy  norm  is  the  one  customarily  used  for  hyperbolic  equations.  In  the  Euler-Bernoulli  equation, 
a  second  derivative  of  in  is  used  to  compute  the  strain  energy.  Such  a  norm  is  used  later  in  the  control 
problems 
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Figure  5:  The  First  Period  of  Trajectory  of  hinged-hinged  beam,  32  Element  Model 


Figure  6:  The  Last  Period  of  Trajectory  of  hinged-hinged  beam,  32  Element  Model 
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DIRK  method  together  have  very  good  accuracy,  but  a  high  mode  beam  needs  high  order 
model.  For  the  first  three  modes,  the  16  element  model  has  given  reasonably  good  accuracy, 
of  order  10“ 4  to  10“*,  but  for  mode  4  and  mode  5,  high  order  mode  is  needed.  For  mode  5, 
the  even  64  element  model  gives  higher  error  than  in  the  lower  mode  case  with  lower  order 
model.  Also,  the  accuracy  slightly  decreases  with  time. 


N 

l|e(<i)Ho 

Mh)\\o 

II«(*3)Ho 

IWMIlo 

Mode  1 

(h=0.01s) 

4 

3.2138c-64 

1.1000c"03 

2.0000c"03 

2.6000c"03 

8 

2.0349c"05 

7.0090c"05 

1.2444c"04 

1.6356c"04 

16 

1.2731c"06 

4.3852c"06 

7.7800c"06 

1.0227c"05 

32 

7.6794c"08 

2.6516c"07 

4.6339e"07 

6.1083c"07 

64 

2.3292c"09 

8.6500c"09 

8.4527c"09 

1.2702c"08 

Mode  2 

(h=0. 0025s) 

4 

6.5000c"83 

1.9300c"02 

8.9200e"°2 

8. 5700 e"82 

8 

4.3288c"04 

1.5000c"03 

5.4000c"03 

6.5000c"03 

16 

2.7080e"os 

9.7218c"06 

3.4042c"04 

4.1599c"04 

32 

1.7090c"06 

6.0937c"06 

2.1271c"05 

2.6029e"05 

64 

9.9915c"08 

3.5991e"07 

1.2721c"06 

1.5493c"06 

(h=0.001s) 

4 

4.6700e~82 

8.5320e"01 

8.5290e"01 

8 

2.1000c"03 

isisilss 

3.6000c"03 

4.7000e_O3 

16 

1.3057c"04 

7.7637c"05 

1.6017c"04 

2.7291c"04 

32 

8.1749c"06 

4.8771c"06 

1.1343c"05 

1.9334c"05 

64 

5.1160e"07 

7.0898e'°7 

1.2094c"06 

Mode  4 

(h=0. 0005s) 

4 

1.2059c"80 

su 

l.lllOe"01 

1.4220e"81 

8 

2.8100e~02 

2.9300e"°2 

3.S100e"01 

16 

1. 7000c"03 

1. 2400c ~02 

1.2000e"°2 

32 

1.0722c"04 

8.2854c"04 

6.9332c"04 

64 

6.7082c"06 

J32BSO 

5.2067e~°5 

4.3146c"04 

Mode  5 

(h=0. 0004s) 

4 

9.0730e"01 

2.3870c"01 

1.0454c"66 

1.7580c"81 

8 

3.8900e"°2 

8.0530e"01 

4.1850e"01 

1.3291c"00 

16 

5.9000c"03 

7.6300e"02 

4.4300e'°2 

4.9700e"°2 

32 

3.6936c"04 

4.9000c"03 

3.4000c"03 

2.2000e_O3 

64 

..  _ 

3.0805c"04 

3.0805c"04 

2.1733c"04 

1.3448c"04 

Table  7:  Lj  norms  of  errors  between  exact  and  numerical  solutions  at  t  =  1,5,8, 10  seconds. 
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Figure  7:  Log  of  Li  error  vs.  log  of  element  number  (t=ls) 


log  o)  Ume 


Figure  8:  Log  of  Li  error  vs.  log  of  time,  with  32  element  model 
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To  obtain  the  energy  norm  (in  Eq.(80)),  the  dwtx  fdx  and  dwnu/dx  need  be  calculated 
first.  The  dwex/dx  can  be  derived  directly  from  exact  solution  Eq.(8l) 

duifT 

— •  =  fico30xcoswt  (82) 

The  numerical  solution  dwnu/dx  can  be  obtained  through  finite-element  method.  First, 
the  local  form  dw^u/dx  is  derived  from  Eq.(7) 

(83) 

ax  i= 1 

Then  the  global  solution  is  assembled  according  to  boundary  condition.  The  relative 
energy  norm  errors  for  the  first  five  modes  are  shown  in  Table  8.  The  logarithms  of  the 
energy  norm  errors  versus  the  logarithm  of  element  number  and  time  are  shown  in  Fig. (9) 
and  Fig.(10). 


Figure  9:  Log  of  energy  norm  error  vs.  log  of  element  number  (t=ls) 


N 

ll«(*i)llo 

NMIIo 

lle(«3)||o 

lk(^)llo  | 

Mode  1 

(h=0.01s) 

4 

3.4000c"03 

3.4000c-03 

1.8000c-02 

1.9200c-05 

8 

3.775 6er°* 

3.6796c-04 

1.2000c-03 

1. 3000c -03 

16 

4.4161e-°5 

4.4161c-05 

8.2840c"05 

8.6877c"05 

32 

5.4688c-06 

5.4596c-06 

6.8722c-06 

7.0547c'06 

64 

6.8015c-07 

6.8039c"07 

6.8710e-07 

6.8496c"07 

Mode  2 

(h=0. 0025s) 

4 

2.3200e-°2 

3.5700c"02 

1.6140e-01 

1.6470c"01 

8 

2.8000c-03 

3.6000c"03 

l.OlOOe'02 

1.2700e-°2 

16 

3.7063c-04 

3.9911c"04 

7.0241c"04 

8.8320c"04 

32 

4.3573e-°5 

4.4436e"05 

5.7915e"05 

6.5879e"05 

j  64 

5.4802c-06 

5.4802c-06 

5.9077c-06 

6.2033c"06 

Mode  3 

(h=0.001s) 

4 

8.4500e-°2 

5.7460e"01 

1.2081e-°° 

1.5980e-°° 

8 

9.3000c-03 

9.3000c-“ 

1. 0000c "°2 

1.0900e"°2 

16 

1.2000c-03 

1.2000c-03 

1.2000c-03 

1.2000e-03 

32 

1.4665c-04 

1.4636c"04 

1.4709c"04 

1.4876c"04 

64 

1.8347c-05 

1.8347c"04 

1.8369e"05 

1.8422e"05 

Mode  4 

(h=0.0005s) 

4 

1.8063c-°° 

2.1125c-00 

1.9710e-01 

2.2720e"01 

8 

4.5800e-°2 

3.2912e-°° 

4.9900e"02 

5.7210e-01 

16 

3.7000c-03 

2.3130c-00 

1.9500e-°2 

1.7700e"°2 

32 

3.9366c-04 

1.4700e"°2 

1.3000c"03 

l.lOOOe"03 

64 

4.4585e-os 

9.2051c'04 

9.1739e-05 

7.6602e~os 

Mode  5 

(h=0. 0004s) 

4 

9.2610e-01 

5.7110c-00 

1.4600c-°° 

7.6780e-°l 

8 

1.7640c-01 

8.8090e-°° 

1.5786e-°° 

1.9580e-°° 

16 

1.0200e-°2 

1.0891e-°° 

7.0200e"°2 

7.2100e"°2 

32 

8.6750c-04 

7.0000e"°2 

5.5000c"03 

3.3000c-03 

64 

9.1475c-04 

4.4000e-03 

3.5354c"04 

2.1217c"04 

Table  8:  Energy  norms  of  errors  between  exact  and  numerical  solutions  at  t  =  1,5,8, 10 

seconds. 


5.4.2  Clamped-free  Beam 

We  consider  the  same  beam  used  in  section  3.2.1  here,  but  with  boundary  condition  fixed  at 
left  end  and  free  at  right  end.  The  mode  shape  function  for  clamped-free  beam  has  following 
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Figure  10:  Log  of  energy  norm  error  vs.  log  of  time  (64  element) 


form 


W»  =  Cn[coshpnx  -  cospnx  -  Csinp\^sinh^nX  ~  sinP"x^  n  =  0, 1, 2, . 

where  Pnl  =  n7r  +  f  +  o(~).  The  values  of  Pnl  for  the  first  five  modes  are 


(84) 


Pxl  =  1.8744 
P2l  =  4.6948 
Pzl  =  7.8540 
P<\  =  10.9956 
Pit  =  14.1372 


The  period  and  energy  preservation  of  DIRK  method  in  this  case  is  examined  in  the  same 
way  as  in  hinged-hinged  beam  case.  The  simulations  have  been  done  on  4,  8,  16,  and  32 
element  models,  all  for  the  first  mode.  The  time  interval  for  each  simulation  is  14,040  seconds, 
which  includes  1111  periods.  The  distances  between  the  first  zero  and  the  2221th  zero  were 
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calculated  based  on  both  analytical  solution  and  numerical  simulation.  The  comparison  of 
the  distances  are  shown  in  Table  9,  where  T  and  T  represent  the  distances  calculated  based 
on  analytical  solution  and  numerical  simulation  respectively. 


N 

T  (sec) 

TXsec) 

T  -  T{ sec)  | 

4 

8 

16 

32 

1.402607614640182C+04 

1.402607614640182e+°< 

1.402610352676545e+M 

1 .40261052296 1071 e+04 

1.402530249457789e+tM 
1.40257401604 2990e+O4 

1 .402575979053479e+o< 
1.402576151727471e+O4 

7. 7365 1 8239325960e "01  I 
3.359859719221276c-01  f 
3.437362306576688c-01  | 
3.437123359999532c-01  | 

Table  9:  Period  Preservation  of  DIRK  method  for  clamped-free  beam 


The  energy  preservation  of  DIRK  method  is  examined  by  comparing  the  norms  of  all 
states  at  initial  time  and  at  the  end  of  the  1111th  period.  The  results  are  shown  in  Table 
50,  which  are  for  4,  8,  16,  and  32  element  models.  In  Table  50,  norm  1  and  norm  2  represent 
the  norms  of  states  at  the  initial  time  and  the  end  time. 


N 

Norm  1 

Norm  2 

Ratio  (n2/nj) 

4 

8 

16 

32 

5.345123709927852e-02 

7.143825321630370e-02 

9.803811033736426c-°2 

1.365134397632705e_O1 

5.344597898707650e"°2 
7.143131244537075e-°2 
9.80336398554456 1 e~02 
1.365072728087120e-01 

9.999016278670547e-01 

9.999028423761716e-01 

9.999544005703163c-01 

9.999548252936177c-01 

Table  50:  Energy  Preservation  of  DIRK  method  for  clamped-free  beam 


5.5  Stabilization  of  a  Beam  by  using  State  and  Output  Feedback 

A  beam  with  a  nonzero  initial  state  (e.g.  nonzero  initial  deflection,  and  /or  initial  deflection 
velocity)  can  be  brought  to  rest  ("stabilized”)  by  using  several  feedback  mechanisms. 

A  clamped-free  beam  can  be  stabilized  by  using  shear  force  control  at  the  tip  which  is 
proportional  to  the  minus  velocity  of  the  tip 

Shear  force  control=  —  &Jju>(f,  /),  t  >  0 
where  k  >  0. 

Alternatively,  the  control  can  be  computed  by  using  state  feedback  found  by  solving  an 
algebraic  Riccati  equation.  The  determination  of  the  feedback  gains  is  possible  only  through 
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approximation.  In  our  setting,  assuming  that  the  finite  element  model  with  N  elements 
provides  acceptable  approximation  of  the  system  state,  the  approximate  Riccati  feedback 
can  be  computed  by  solving  an  algebraic  Riccati  equation  for  the  finite  dimensional  model. 

For  problems  with  unbounded  control  operators,  the  question  of  convergence  of  such  a 
procedure  when  JV  — ►  oo  is,  in  general,  covered  by  Lasiecka  and  Triggiani  [36]  The  case  of  an 
Euler- Bernoulli  beam  is  in  the  context  of  Linear- quadratic  regulator  has  also  been  discussed 
by  Gibson  and  Adamian  [22],  however  their  setup  involves  an  attachment  of  the  beam  to  a 
hub  at  x  =  0  which  makes  the  input  operator  bounded. 

The  convergence  conditions  given  by  Lasiecka  and  Triggiani  require  that  certain  inequal¬ 
ities  be  satisfied  by  the  approximation  method  uniformly  for  all  N  -*  oo.  Whether  these 
inequalities  hold  in  the  case  of  a  cubic  finite  element  method  appears  to  be  an  open  prob¬ 
lem.  Even  if  they  hold,  pointwise  converge  of  distributed  gain  feedbacks  is  not  guaranteed. 
However,  even  without  the  theoretical  proof  of  convergence,  the  cubic  finite  element  method 
appears  to  produce  numerically  convergent  state  trajectories  as  N  -*  oo. 

The  cost  function  used  in  simulations  is  of  the  form 

JQ  (II  llff*  +  II  Wi^dt  +  R  J  u2(t)dt  (85) 

or 

A00 (II  w(*.  •)  Hi2  +  II  •)  111 »)dt  +  rI°°  «a(0*  (86) 

JO  Jo 

In  either  case,  the  approximation  by  finite  elements  reduces  these  cost  functions  to  fa¬ 
miliar  quadratic  cost  functionals  of  the  finite-dimensional  control  theory 

J  =  Jq  (xT(t)Qx(t)  +  uT(t)Ru(t))dt  (87) 

where  x(t)  =state  vector,  u(t)  =control  vector.  In  the  beam  case  with  energy  norm,  the 
matrices  Q,R  are 


L 

0  M  I’ 


R  =  r 


(88) 
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In  actual  computations,  the  mass  matrix  Af  is  split  by  Cholesky  decomposition 


M  =  McMj 

and  the  state  vector  is  transformed  accordingly 


x  = 


r  mt  o 

0  Mj  , 


S/m 


Hence  the  matrix  Qm  actually  used  in  computation  had  the  form 


(89) 


(90) 


where 


Qm  = 


Lm  O' 

0  I  ’ 


R  =  r(scalar) 


Lm  =  M~1L(Me)~1 


(91) 

(92) 


In  the  case  of  the  second  cost  function,  the  matrix  Qm  is  an  identity  matrix,  because  it 
arises  from  a  Cholesky  transformation  applied  to  Q  =  diag{M,M}. 


5.5.1  Clamped-free  Beam 

Parameters  used  in  simulation  were  /  =  3m,  El/p  =  0.1.  The  initial  conditions  were 
assumed  in  a  form  of  modal  shape  functions 


wn(x)  =  cosh(0nx)  -  cos(0nx)  -  Ri(sinh(/3nx)  -  sin(f)nx))  (93) 


where  Rx  =  {cosh(Pnl)  +  cos(/3nl))/(sinh(0nl)  +  sin(0nl)) 

and  /?„/  =  [1.8744,4.6948,7.8540,10.9956,14.1372]  for  (n=l,...5).  Parameter  pn  relates 
to  modal  frequencies  by  the  formula 


(94) 


Figures  (28)-(31)  show  the  distribution  of  eigenvalues  corresponding  to  R  —  100,  N  = 
8, 16, 32  and  N  =  100.  Characteristically,  most  of  the  eigenvalues  are  distributed  along  a 
vertical  strip  parallel  to  the  imaginary  axis.  This  suggests  that  the  LQR  control  through 
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shear  force  has  an  effect  analogous  to  viscous  damping.  The  bending  of  the  eigenvalue  locus 
towards  the  imaginary  axis  at  high  frequencies  occurs  at  the  frequencies  which  are  artifacts 
of  the  approximation.  One  can  see,  by  comparing  eigenvalue  plots  for  N  =  8, 16, 32, 100  that 
as  N  increases,  the  number  of  eigenvalues  which  remain  in  the  vertical  strip  increases.  Thus, 
the  eigenvalues  which  for  low  N  were  seen  approaching  the  imaginary  axis  are,  for  high  N, 
seen  to  be  distant  from  the  imaginary  axis.  It  is,  therefore,  likely,  that  as  N  — ►  oo  all  the 
eigenvalues  of  large  modulus  still  be  in  the  vertical  strip  bounded  away  from  the  imaginary 
axis.  This  would  suggests  the  absence  of  a  "spillover”  in  case  of  a  Riccati  feedback. 

The  distance  from  the  set  of  eigenvalues  to  the  imaginary  axis  depends,  for  a  fixed  N,  on 
the  parameter  R  in  the  cost.  Figures  (32), (33)  show  that  distance  as  a  function  of  parameter 
R,  for  N  —  8  elements.  It  is  interesting  to  observe  the  existence  of  a  sharp  minimum  of 
max(real(A,)),  that  is  the  existence  of  a  maximum  of  the  distance  functions.  This  occurs  for 
Rm  =  1.64  and  yields  the  maximum  decay  rate  Amar  =  0.87. 

Although  one  would  be  tempted  to  think  that  /£*  yields  the  "best”  transients  to  zero, 
this  is  not  necessarily  so.  In  fact,  by  comparing  Figures  (22)-(24)  one  can  see  that  the  best 
transients  occur  for  R  approximately  equal  100. 

Figure  (25)  shows  a  three  dimensional  view  of  the  beam  stabilization  process  with  the 
initial  function  being  the  first  modal  shape. 

For  output  feedback,  Lagnese  [33]  provides  a  formula  for  the  decay  rate  in  an  upper 
bound  on  the  energy  norm.  This  formula,  after  a  time  scale  conversion,  yields  the  best  decay 
rate 


.  i  [eI 

"  2L»fl  +' W)  VT  (95) 

For  the  specific  beam  parameters,  this  yields 

Amar  =  0.01029  (96) 

This  appears  to  be  a  very  conservative  estimate.  Simulations  show  that  much  better 
decay  is  possible  with  output  feedback. 
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Figure  (26)  give  comparison  of  transients  under  near  optimal  (numerically)  value  of  the 
feedback  gain  k  =  0.022  in  the  output  feedback,  and  the  state  feedback  corresponding  to 
R  =  100.  The  state  feedback  is  seen  to  provide  about  twice  as  fast  stabilization  as  the 
output  feedback.  This  suggests,  that  an  improvement  in  the  output  feedback’s  performance 
is  possible  if  one  adds  sensors  to  the  beam,  so  as  to  provide  dynamically  the  information 
about  the  beam’s  shape.  Alternatively,  one  can  insert  an  observer  in  the  loop. 

5.5.2  Hinged-hinged  Beam 

While  a  hinged-hinged  beam  with  a  torque  control  at  one  end  is  less  likely  to  be  encountered 
in  typical  applications,  it  represents  an  interesting  case  mathematically.  Without  control, 
it  is  a  problem  with  symmetry  induced  by  the  boundary  conditions,  thus  offering  easier 
interpretations  that  the  unsymmetric  clamped-free  and  hinged-free  problems. 

The  beam  investigated  numerically  has  parameters  /  =  1,  El  j p  —  0.02.  Figures  (40)- 
(44)  show  the  closed-loop  eigenvalue  patterns  for  the  system  under  state  (Riccati)  feedback 
for  A  =  8, 16, 32. 

In  contrast  to  the  clamped-free  beam,  the  closed-loop  eigenvalues  are  now  distributed 
primarily  within  the  conical  sector  in  the  left  half  plane.  The  eigenvalues  of  large  modulus 
are  larger  negative  real  parts.  This  pattern  is  consistent  for  all  values  of  N,  and  R,  and  is 
distributed  only  by  the  few  eigenvalues  at  the  end  of  the  chain.  For  very  low  values  of  R  one 
also  finds  two  large  real  eigenvalues  for  in  the  left  half-plane. 

The  transient  corresponding  to  the  state  feedback  is  shown  on  Figure  (27). 

5.6  Repositioning  of  a  Beam  using  State  or  Output  Feedback 

A  beam  can  be  repositioned  by  using  a  external  input,  torque  or  shear  force,  at  either  end 
of  beam,  and  stabilized  via  output  feedback  or  state  feedback  compensator.  In  the  latter 
case,  the  states  of  system  are  assumed  to  be  known,  and  state  estimation  will  be  discussed 
in  next  section.  An  external  input,  a  torque  as  bending  moment  or  a  force  as  shear  force,  is 
applied  at  either  ends,  depending  on  the  boundary  conditions  of  beam.  Only  is  static  input 
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considered  in  the  simulation,  with  a  proportional  gain  to  drive  a  beam  to  desired  position. 
The  output  feedback  compensator  is  based  on  feedbacks  in  velocity-type  term  and  position 
or  slope.  The  state  feedback  compensator  is  designed  by  using  operator  Riccati  equations. 
In  addition,  a  Posicast  type  pre-compensator,  or  preshaping-command  input  technique,  is 
used  to  reduce  the  vibration  of  beam,  due  to  the  forcing  function. 

In  this  simulation,  FE  model  is  used  to  approximate  the  PDE  model  of  Euler-Bernoulli 
beam,  and  the  second-order  DIRK  method  is  used  as  integrating  method.  Undamped  hinged- 
free,  and  clamped-free  beams  are  considered  here. 


5.6.1  Input  Gain 

As  mentioned  before,  the  external  inputs  for  beams  in  the  simulation  are  a  torque  as  bending 
moment  applied  at  left  end  for  hinged-hinged  beam,  and  a  shear  force  applied  at  right  end 
for  clamped-free  beam.  In  either  case,  the  external  input  r  is  a  step  function  with  amplitude 
G,,  i.e.  r  =  Gjl(t),  which  is  adjusted  to  the  desired  value  of  the  system  output.  The  G,- 
can  be  obtained  through  the  static  equation  of  beam  and  desired  values.  The  Gi  can  be 
interpreted  as  a  gain  needed  at  the  system  input  to  produce  the  desired  output  yj. 

First,  the  input  gain  G,-  for  output  feedback  design  is  considered.  The  undamped  system 
with  output  feedback  can  be  described  as 


y 


u(<) 


-Lm  y  +  Qmu(i) 


-KC 


+  r 

-'**i [ co  ft] [5] 

—K\C\y  —  K2C2y  +  r 


+  r 


(97) 

(98) 


where  K\,  C\  and  K2,  C2  are  feedback  gain  vectors  and  output  matrices,  corresponding 
to  the  variables  y  and  y  respectively.  Substituting  Eq.(97)  into  Eq.(98)  and  rearranging  it 
leads  to 


y  =  ~(Lm  +  QmmKiCx)y  -  QmK2C2y  +  Qmr 


(99) 
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The  static  equation  is  obtained  by  setting,  for  t  —  oo,  y  =  0,  y  =  0,  and  y  =  yd  in  Eq.(99) 


as 


0  =  ~(Lm  +  QmmKxCx)yd  +  QmG4(<)  (100) 

where  yd  is  the  desired  value.  Hence,  the  value  of  G,-  can  be  obtained  by  solving  Eq.(lOO). 


G;  = 


QmiLm  +  QmK,Cx)yd 


QmQm 

Similarly,  for  state  feedback  the  value  of  G,-  can  be  calculated  by 


(101) 


Gi  = 


Qm(Lm  +  QmKi)yd 

QlQm 


(102) 


where  Kx  is  state  feedback  gain  corresponding  to  the  variable  y. 


5.6.2  Output  Feedback  Control 

The  simulations  for  output  feedback  control  were  done  in  following  cases: 

(i)  hinged-free  beam  with  a  control  torque  at  left  end,  and  feedback  on  the  velocity  of 
slope  at  left  end,  plus  slope  at  left  end,  or  position  at  right  end.  However,  whether  the  latter 
two  feedbacks  are  needed  depends  on  the  flexibility  of  beam. 

(ii)  clamped-free  beams  with  a  shear  force  input  at  free  end,  and  a  feedback  as  the  velocity 
of  position  at  free  end  (the  "robotic  arm”  configuration). 

In  the  simulation,  system  time  responses,  corresponding  to  the  flexibility  of  beam  and 
element  number  of  FE  models  was  studied.  Also  the  comparison  of  output  feedback  stabi¬ 
lization  between  hinged-free  and  clamped-free  beam  was  done  in  the  term  of  feedback  gain, 
and  stability  margin.  A  lightly  damped  beam  was  used  in  the  simulation  with  parameters 
as  /  =  3,/?  =  0.1,  k  =  0.02 ,EI  =  0.02  to  100.  The  desired  static  position  at  the  tip  of  beam 
is  yd  =  1. 

The  time  responses  of  position  at  tip  for  a  hinged-free  beam  with  flexibility  El  = 
0.1, 0.5, 1  and  100  are  shown  in  Fig.(ll).  In  this  case,  three  sensors  were  assumed,  mear 
suring  the  slope  at  left  end,  position  at  right  end,  and  velocity  of  slope  at  left  end,  with 
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feedback  gain  k  =  [—0.50.05  —  1.5].  The  beam  waa  repositioned  from  zero  to  one  at  tip.  A 
dynamic  trajectory  for  repositioning  of  a  hinged-free  beam,  displacement  vs  time  coordinates, 
is  shown  in  Fig.(12). 


Figure  11:  Repositioning  of  hinged- free  beams  with  output  feedback 

From  Fig.(ll  ),  we  can  see  that  the  beam  with  smaller  El  has  longer  settling  time  and 
bigger  overshot.  The  beam  with  El  =  100  reach  static  state  in  about  7  seconds,  while 
the  beam  with  El  =  0.1  could  not  settle  down  at  all  during  20  second  simulation  interval. 
The  feedback  gains  used  in  the  four  cases  were  the  same.  Adjusting  the  feedback  gains 
according  to  different  beams  would  improve  dynamical  trajectory  somewhat.  At  least  one 
more  feedback  besides  the  velocity  of  slope  at  left  end,  is  needed  to  stabilize  the  system. 
The  reason  for  this  is  that  hinged-free  beam  has  a  double  eigenvalue  at  zero,  one  of  which 
cannot  be  moved  to  the  left  half  plane  by  just  one  feedback. 
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Figure  13:  Repositioning  of  a  hinged-free  beam  with  a),  output  feedback  of  tip  position, 
slope  and  slope  rate,  b).  state  feedback  with  observer.  Much  better  decay  of  vibration  is 
achievable  in  case  b. 

dynamically  estimated  from  available  sensor  measurements.  This  task  will  be  performed  by  a 
dynamic  observer  -  i.e.  a  state  estimator  for  the  structure,  see  next  section.  The  dynamically 
changing  reconstructed  state  is  then  fed  into  the  state  feedback  controller. 

The  state  feedback  designed  by  discrete  approximations  and  LQG  methodology  leads  to 
an  interesting  pattern  of  closed  loop  eigenvalues  (Fig.(14)).  Namely,  the  optimally  placed 
closed  loop  eigenvalues  form  a  conical  pattern  along  straight  lines  going  into  the  left-hand 
plane.  The  high  frequency  modes  exhibit  much  larger  damping  than  the  low  frequency 
modes.  This  conclusion,  already  seen  in  the  case  of  a  hinged-hinged  beam,  is  thus  also  true 
for  a  hinged-free  beam. 
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Figure  14:  A  typical  example  of  closed  loop  eigenvalues  with  Riccati  feedback  design 
5.6.4  Feedback  Control  with  Posicast  Pre-compensator 

Control  systems  built  by  using  output  feedback  respond  poorly  to  external  perturbations  and 
have  modest  stability  margins.  In  addition,  such  systems  poorly  execute  slewing  command 
maneuvers  (see  Fig.(12)  and  (13)).  While  the  latter  can  be  improved  by  command  signal 
shaping  using  low-pass  prefilters,  Posicast  pre-compensators,  or  maneuver  pre-programming, 
the  weak  damping  and  small  stability  margin  remains  the  basic  feature  of  such  control 
systems. 

Posicast  is  one  of  preshaping  command  input  techniques  which  are  used  to  reduce  system 
vibration.  It  involves  breaking  a  step  of  a  certain  magnitude  into  two  smaller  steps,  one  of 
which  is  delayed  in  time.  This  results  in  a  response  with  a  reduced  settling  time.  In  effect, 
superposition  of  the  responses  leads  to  vibration  cancellation.  Posicast  technique  used  here 
is  designed  based  on  the  dynamic  characteristics  of  the  closed-loop  system,  which  is  shown 
in  Fig.(15). 
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Figure  15:  Preshaping  the  command  input  with  Posicast  Filter:  r(t)-command  to  position 
the  tip  of  the  beam;  F-feedback  (state  or  output) 


In  the  beam  problem,  Posicast  pre-compensator  is  used  to  preshape  the  command  so 
that  the  least  damped  modes  in  a  beam  won’t  appear  in  the  response.  In  order  to  do  this, 
P(s)  should  have  zeros  at  the  least  damped  poles  of  the  closed-loop  system.  Assume  that 
the  eigenvalues  of  closed-loop  system  at  the  least  damped  poles  are 


An  —  £n  i  n  —  1,2, ...,m 


(103) 


The  Posicast  pre-compensator  has  its  transfer  function  corresponding  to  each  individual 
mode  as 


Wn(s)  = 


1 


(104) 


1  +Kn  1  +Kn 

where  Kn  is  fractional  overshoot  of  the  step  response  for  the  system  to  be  controlled  and 
Tn  is  the  time  to  the  peak  for  the  step  response.  We  can  express  Kn  and  Tn  as 


(105) 

(106) 


The  transfer  function  has  zeros  at  \n  in  Eq.(103).  Cascaded  Posicast  compensator  in 
multi-mode  case  has  a  form 
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P{*)  =  ft  Wk(s)  (107) 

k=l 

The  Posicast  pre-compensator  has  zeros  at  A*,  k  =  1,2,  ...m.  The  repositioning  of  a 
hinged- free  beam  with  output  feedback  using  Posicast  pre-compensator  is  shown  in  Fig.(16), 
together  with  the  response  not  using  Posicast.  The  Posicast  pre-compensator  used  here  was 
built  based  on  the  first  three  modes. 
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Figure  16:  Repositioning  of  a  hinged-free  beam  (El  =  0.5)  through  output  feedback  with 
and  without  Posicast  pre-compensator 

The  simulation  shows  that  the  oscillation  of  repositioning  has  been  eliminated  and  settling 
time  is  quite  short  in  this  case. 

5.7  Observer  Design 

One  of  the  main  steps  in  the  construction  of  the  feedback  operator  for  the  boundary  control 
of  the  beam  is  the  construction  of  the  observer  which  will  compute  the  state  of  the  beam 
(deflection  and  it’s  time  derivative)  at  every  time  instant  t  based  just  on  the  measurement  of 
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one  scalar  valued  function  z(t).  Regular  Luenberger  observer  in  discrete-time  and  adaptive 
observer  designs  are  discussed  in  this  section. 

5.7.1  Regular  Observer  in  Discrete-time 

The  main  numerical  difficulty  here  is  that  the  observer  has  an  internal  feedback  loop  which 
shifts  the  eigenvalues  of  the  A  matrix  by  transforming  A  into  A  —  LC ,  where  C  is  the  output 
matrix,  and  L  is  the  observer  gain.  The  matrix  A  —  LC  corresponds  to  system  that  is 
no  longer  purely  oscillatory,  and  thus  does  not  lend  itself  to  the  use  of  the  DIRK  method 
for  second  order  ODE’s.  One  then  has  to  use  a  method  suitable  for  damped  oscillatory 
problems,  which  puts  us  back  into  the  consideration  of  the  first  order  state  equations  system 
which,  as  discussed  before,  suffers  from  conditioning  problems.  We  are  interested  in  having 
an  observer  of  reasonably  low  dimension,  yet  capable  of  reproducing  the  motions  of  several 
principal  modes  of  the  beam. 

In  the  process  of  designing  the  observer  (i.e.  finding  the  gain  matrix  L)  one  has  a  choice 
of  first  designing  L  for  a  continuous  time  model,  and  then  applying  the  discretization  in  time 
to  matrix  A  —  LC,  or  alternatively,  one  can  first  compute  the  discrete-time  equivalent  of 
matrices  A ,  C ,  and  then  perform  the  observer  design  in  the  discrete-time  framework.  We 
found  that  the  two  steps  do  not  commute,  and  that  the  second  choice  is  a  more  accurate 
approach.  The  first  choice  requires  using  a  predicted  value  of  z(t  +  hi),  where  hi  is  a  fraction 
of  step  size  h  dictated  by  the  DIRK  method. 

The  FE  model  of  a  beam  in  discrete-time,  after  using  DIRK  method,  can  be  described 
as 


a?n+i  =  Fxn  +  Gun  (108) 

zn  =  Hnixn  (109) 

Designing  a  observer  based  on  the  discrete-time  model,  we  have 

*n+l  =  Fxn  +  L(zn-  zn)  (110) 

Zn  =  Hn'Xn  (111) 
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where  rii  and  n2  can  be  different  in  the  case  that  observer  has  lower  dimension  model. 
The  observer  gain  L  can  be  obtained  by  shifting  eigenvalues  left  or  using  Riccati  equation. 

The  observer  errors  for  a  hinged-hinged  beam  in  the  first  two  modes  are  shown  in  Fig.(17 
)  and  Fig.(18  ).  The  beam  was  approximated  by  64  element  FE  model,  and  the  observers 
were  in  lower  dimension,  8  element.  The  observer  was  built  based  on  only  one  measurement, 
the  slope  of  beam  at  right  end.  The  computational  experiments  showed  that  an  entire  shape 
of  the  beam  can  be  rapidly  reconstructed  while  the  beam  vibrates,  with  good  accuracy  for 
the  first  few  natural  modes.  The  number  of  modes  accurately  reconstructed  can  be  increased 
at  the  expense  of  computational  effort. 

5.7.2  Adaptive  Observer  Design 

The  adaptive  observer  for  Euler- Bernoulli  beam  is  designed  to  estimate  the  states  and  the 
parameters  of  system  in  FE  model,  with  unknown  initial  conditions,  through  only  a  few  mea¬ 
surements.  The  unknown  parameters  can  be  only  the  coefficients  of  stiffness  and  damping. 
Considering  the  unknown  parameters,  the  FE  model  for  beam  can  be  expressed  as 

y  =  -qiKmy  -  q2Lm  +  Qm  (112) 

where  qi  —  k/p ,  q2  =  El  f  p  are  unkown  parair  >ers,  and  Km  =  M~lKq,  Lm  =  A/-% 
and  Qm  =  M;XQ  with  Mq  =  M/q,  Kq  =  K/k  and  Lq  =  L/EI  (M,K,L  see  Eq.(17)). 
Converting  the  equation  (112)  into  a  standard  state-space  form,  we  have 


x 

y 


0 

I 

i  + 

0 

—q2Lm 

-q\Km 

Qm  . 

(113) 

(114) 


However,  the  typical  adaptive  observer  design  schemes  appeared  in  literatures  are  based 
on  some  specific  system  representations,  for  example,  the  system  described  by  a  differential 
equation  of  the  form 
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(115) 


x  =  [a|A}z  +  bu 
y  =  cTx  =  xi 


where  a  =  [aijaj;  ...;an]r,h  =  f&i;  63;  ...&n]T,  the  (n  x  (n  —  l))matrix  A  is  known  and 
(F  =  [1;0;  ...0].  The  vectors  a  and  b  represent  the  unknown  parameters  of  the  system.  The 
matrix  A  is  in  the  canonical  form 


(116) 


where  I  €  72.(n-1)x(w“l)  is  an  identity  matrix,  and  A,-  are  the  eigenvalues  of  system. 

Such  adaptive  observer  design  scheme  is  not  suitable  for  our  beam  problem  because  of 
the  following  reasons 

•  The  transformation  from  FE  model  to  specifical  form  causes  poor  numerical  condition, 
specially  for  large  dimension  model. 

•  State  variables  in  the  specifical  representation  no  longer  has  original  physical  interpre¬ 
tation,  and  recovery  original  state  variables  through  transformation  is  too  complicated. 

•  There  are  2 n  parameters  to  be  tracked  even  though  only  2  parameters  in  FE  model  are 
unknown. 

An  adaptive  observer  design  scheme  based  on  state-space  FE  model  of  E-B  beam  is 
considered  here,  which  using  sensitivity  function  concept  to  define  adaptive  updating  law. 
For  the  system  described  in  Eq.(113),  the  observer  can  be  designed  as 


x  = 
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-qjLm 


I 

-qiK„ 


x  +  L(y-y)  + 
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(118) 
(119) 
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where  qx  =  k/p  and  92  =  El /p  are  unknown  parameters.  To  obtain  updating  law  for 
the  observer,  the  sensitivity  functions  of  the  system  with  respect  to  parameters  91  and  93 
are  studied.  Sensitivity  functions  provide  first-order  estimates  of  the  effect  of  parameter 
variations  on  solutions.  The  observer  equation,  Eq.(117),  can  be  written  to 


x  =  A(q)x  +  L(y  -  y)  4-  Bu  (120) 

where  9  =  [91  92]^,  and  let  input  u  be  zero. 

Define 

«•)  =  ^  (121) 

The  sensitivity  function  around  nominal  solution  9  is  expressed  as 

{(l)  =  ^(«X(0  +  ^i(f.»)  (122) 

We  have 


6  = 
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-q-iLm 
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£1  + 

'  0  o' 

.0  -Km 

x(t) 

(123) 
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(124) 

The  adaptive  updating  law  can  be  described  as 


9i 

=  — 7eiC£i 

(125) 

92 

=  -7txCii 

(126) 

9i(0) 

=  9*x 

(127) 

92(0) 

=  9.2 

(128) 

where  ei  =  y  —  y  is  observer  error,  and  7  is  speed  coefficient. 
The  numerical  simulation  procedure  for  observer  is  as  following 
Step  1:  a)  Given  initial  parameters  9,-,,  9;,,  and  measurement  y. 
b)  Calculating  x  and  y  by  solving  Eq.(117)  and  (118). 
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Step  2:  Computing  sensitivity  functions  £1  and  by  using  Eq.(123)  and  (124). 

Step  3:  Updating  parameters  41  and  92  by  using  Eq.(125)  and  (125).  With  updated 
parameters,  go  back  to  Step  1-b)  and  repeat  the  whole  procedure  until  the  observer  error  is 
small  enough. 

The  observer  gain  L  in  Eq.(117)  can  be  obtained  by  placing  poles  on  the  system  corre¬ 
sponding  to  each  updated  parameters.  However,  doing  pole  placement  in  each  updating  loop 
brings  calculation  difficulties,  especially  when  there  are  multiple  eigenvalues  in  system.  What 
we  used  in  the  simulation  is  to  design  L  initially  for  initial  parameters  by  pole  placement, 
and  redesign  it  only  when  closed-loop  system  A(q)  —  LC  become  unstable. 

The  simulation  results  for  hinged-hinged  beam  with  parameters  91  =  0.2  and  92  =  1  are 
shown  in  Fig.(20)  and  Fig.(21  ).  The  output  error  converged  quite  good  in  30  seconds  even 
one  of  the  parameter  (91)  did  not  really  converge  to  its  truth  value.  The  Fig.(21  )  shows  the 
decay  of  the  state  error  norm  for  all  states  and  the  states  representing  deflection. 

6  Computational  Experiments  with  a  Timoshenko  Beam 
Model 

In  addition  to  the  experiments  with  the  Euler-Bernoulli  beam,  two  finite-element  models 
were  developed  for  the  Timoshenko  beam  model. 

As  is  well  known,  the  Timoshenko  model  taken  into  account  the  shear  deformation  of 
beam’s  filaments.  This  leads  to  a  single  fourth  order  in  time-fourth  order  in  space  PDE ,  or 
alternatively,  to  a  coupled  system  of  two  second  order  in  space-second  order  in  time  PDEs. 

The  Timoshenko  model  is  known  to  result  in  a  lower  and  more  realistic  rate  of  growth 
of  natural  beam  frequencies  that  exhibited  by  the  Euler-Bernoulli  model.  A  number  of 
experiments  performed  by  other  researchers  in  labs  showed  good  agreement  of  experimen¬ 
tally  determined  frequencies  with  the  computed  ones.  The  Timoshenko  beam  actually  has 
two  chains  of  eigenvalues,  one  corresponding  to  the  deflection  modes,  the  other,  at  higher 
frequencies,  to  the  shear  deformation  modes. 

Our  analysis  to  date  showed  that  the  fourth  order  Timoshenko  model  exhibits  certain 


57 


peculiar  phenomena  in  the  finite  element  context,  namely  the  eigenvalues  of  cubic  FE  approx¬ 
imation  are  not  necessarily  restricted  to  the  imaginary  axis.  This  can  be  shown  to  be  caused 
by  an  inaccurate  reproduction  of  eigenfunctions  by  the  approximating  system,  especially  at 
high  frequencies. 

The  second  order  Timoshenko  system  is  free  from  this  inconvenience,  however  its  use 
with  cubic  elements  entails  twice  as  large  dimension  of  the  approximating  system  as  the 
fourth  order  model. 

We  have  obtained  a  good  approximation  of  natural  frequencies  using  even  a  piece-wise 
linear  system  of  finite  elements  and  a  second  order  Timoshenko  system.  Details  of  this 
research  will  be  described  in  another  report  due  to  that  research  being  done  beyond  the 
termination  date  of  the  present  subcontract. 

Our  current  work  on  this  project  involves  combining  single  beams  into  an  interconnected 
system  of  Timoshenko  beams,  as  described  by  Lagnese  et  al  [34]  This  work  will  enable  us 
to  model  experimental  structures,  such  as  e.g.  the  CSI  Phase  2  structure  at  NASA  Langley. 
This  work  will  be  described  in  subsequent  reports. 

7  Software  Developed  Under  This  Project 

Software  developed  under  this  project  was  done  in  Matlab,  version  3.5  and  version  4.0.  These 
are  programs  corresponding  to  various  beam  equations  (Euler  Bernoulli  and  Timoshenko,  the 
latter  in  both  fourth  order  scalar  version  and  in  the  second  order  system  version),  boundary 
conditions,  type  of  feedback  control,  and  problem  to  be  .solved. 

One  of  the  critical  issues  in  developing  these  programs  was  the  large  dimension  of  the 
system  matrices.  This  required  avoiding  the  standard  Matlab  numerical  integration  routine 
”lsim.m”,  which  easily  bogs  down  on  large  FE  models  (the  apparent  reason  for  this  is  a 
sensitive  Pade  approximation  of  eAt  used  by  lsim.m)  .  However,  replacing  lsim  by  a  custom 
designed  iterative  scheme  could  lead  to  a  very  slow  computation,  because  Matlab  gets  slowed 
down  by  a  large  number  of  iterations  within  loops.  This  difficulty  was  overcome  by  using 
the  DIRK  method  to  produce  matrices  of  a  discrete-time  system  (essentially,  we  use  DIRK 
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only  to  compute  an  energy-conserving  approximant  to  the  exponential  matrix  eAl),  and  then 
using  a  high  speed  Matlab  kernel  "ltitr.m”  to  solve  the  discrete-time  linear  system  without 
using  any  Matlab  loops.  This  turned  out  to  be  a  fast  procedure. 

A  brief  list  of  selected  programs  is  given  below. 

bmhin.m  -  computes  the  finite  element  model  for  a  hinged-hinged  E-B  beam,  along  with 
the  modal  frequencies. 

bmcla.m  -  computes  the  finite  element  model  for  a  clamped- free  E-B  beam,  along  with 
the  modal  frequencies. 

dsl*.m  -  these  programs  integrate  the  state  trajectory  of  the  beam’s  FE  model  by  using 
the  second  order  or  the  first  order  DIRK  algorithm  (*  denotes  the  wild  card  character,  there 
axe  several  version  of  this  program) 

dslctr*.m  -  these  programs  compute  the  trajectory  of  the  closed  loop  system  under  either 
the  Riccati  feedback  or  the  velocity-position  feedback. 

wm.m  -  this  program  computes  the  polynomial  interpolations  between  the  mesh  points, 
allowing  for  a  graphic  three  dimensional  display  of  beam’s  motion  in  time. 

dslobs*.m  -  these  programs  solve  the  observer’s  equations  for  the  E-B  model,  given  the 
placement  of  sensors  on  the  beam  and  the  desired  shift  of  the  eigenvalues  into  the  left  half 
plane. 

adobs*.m  -  these  programs  solve  the  adaptive  observer  problem  for  a  beam  with  unknown 
parameters  El  and  p,  by  computing  the  solutions  of  the  sensitivity  differential  equations, 
and  by  using  a  variation  of  the  gradient  rule  for  the  adaptive  observer  algorithm. 

Currently,  software  under  development  includes  a  FE  model  for  a  rectangular  Kirchoff 
plate  using  cubic-quartic  FEs,  and  a  model  for  a  system  of  interconnected  Timoshenko 
beams  with  deflection  and  axial  deformation  modeled  by  piecewise  linear  FEs  in  the  system 
of  second-order  Timoshenko  equations. 

derivtim*.m  -  these  programs  compute  the  fourth  order  in  time  Timoshenko  FE  model 
using  cubic  Hermite  splines. 
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A  Appendix:  Derivation  of  the  Finite  Element  Model 

A.l  Variational  formulation 

Euler-Bernoulli  beam  is  described  by 


dM‘  x)  tgM(t,x)  E,d^x)  =  0 
Hdt2  dt  dx* 


(129) 


where  0  <  x  <  /,  El  is  the  product  of  the  modulus  of  elasticity  and  the  moment  of 
inertia,  p  =  A-m  and  A  is  the  cross  sectional  area  of  beam  and  m  is  mass  per  unit  volume, 
and  k  is  the  damping  coefficient  of  beam.  w(t,x)  is  the  vertical  deflection  of  beam. 

Before  using  finite  element  approximation,  a  variational  formulation  is  formed  which  re¬ 
casts  a  given  differential  equation  in  an  equivalent  integral  form  by  trading  the  differentiation 
between  a  test  function  and  the  dependent  variable.  We  construct  the  variational  form  of 
Eq.  (129)  over  element  e,  with  test  function  v 


where 
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A. 2  Finite  Element  Method 


We  divide  a  beam  into  N  element  with  N  -f-  1  nodes,  as  shown  in  Fig.(l),  and  assume  that 
displacement  w  over  each  element  is  interpolated  by  an  expression  of  the  form 


where 


w  =  #«!■' + *<■>«<•> + 
j= i 

u$*>  =  =  u>(xe),  =  itj**  ss  to(xe+i)  as  deflections, 

u?  —  ^  ~  0(xe),  —  0^  =  9{xe+i)  as  slopes. 


and  •p<*)  are  interpolation  functions 


4'1  =  + 


A"  =  -(x-x,)(l-^)3 


A’1  =  3(i^)’-2(^)3 
A‘>  =  -(x-X.)((l=i)’-^] 


(132) 


(133) 


The  interpolation  implies  that  at  any  arbitrarily  fixed  time  t  >  0,  the  function  w  can  be 
approximated  by  a  linear  combination  of  with  being  the  value  of  w  or  dw/dx,  at 
time  t,  at  both  end  of  element  e.  In  other  words,  the  time  and  spatial  variations  of  w  are 
separable.  Substituting  for  v  =  4>j'\x)  and  Eq,(132)  into  variational  formulation  Eq.(130), 
we  obtain 


Jfx*+ JL  (Pm  d*<5,- 


■ir  dUj 


(134) 


[Af<e>]{u}  +  [tf(e)]{u}  +  [I(e)]{u}  =  Q(e) 


(135) 
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where 


u  =  ( Ui  ti2  u3  u4  ]T 
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dx2  dx2 

For  the  case  in  which  the  p,  A:  and  El  are  constant  over  an  element,  the  element  matrices 
Af(«)t  k(')  and  £,(*)  are  given  by  (these  are  called  the  mass,  damping,  and  stiffness  matrix, 
respectively) 
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A. 3  Assembly  of  Element  Models 


The  global  model  is  obtained  by  assembling  all  element  models,  in  the  way  that  the  coeffi¬ 
cient  matrices  of  element  models  overlap  at  all  intermediate  nodes.  We  note  the  following 
correspondence  between  the  local  variables  Uj  and  the  global  variables  Uj  (see  Fig.(19)): 
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Figure  19:  Assembly  of  the  beam  element,  N  =  2 
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Since  there  are  two  global  variables  ({/,)  per  node,  the  coefBcients  associated  with  re¬ 
peated  local  variables  should  add  up.  In  other  words,  the  global  coefficients  have  contri¬ 
butions  from  both  consecutive  elements.  For  example,  in  two  element  beam,  the  global 
coefficients  K33,  K34,  K43  and  K44  have  contributions  from  both  element  1  and  2. 


is  _  1  fc'l2)  is  _  fcdri  1  isC*) 

-ft 33  —  ft  33  +  A11  -034  —  ft  34  +  ft  12 

K43  =  Kg  +  Kg  K44  =  Kg  +  Kg 

In  general,  the  assembled  matrices  for  the  assembly  of  beam  element  have  the  form  shown 
in  Eq.(145).  The  global  model  is  in  the  form 

MU  +  KU  +  LU  =  Q  (143) 
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where  coefficient  matrices  A/, K  and  L  are  in  (2 N  +  r)  dimension,  which  fire  assembled 
in  the  way  as  shown  in  Eq.(145),  and  the  boundary  condition  Q,  in  (r  =  0, 1  or  2)4has  the 
form 
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where 


Fo  =  «','>  =  £(*/£)]_  o 

«.  =  < 3S”  =  [£/0).=o 

=  -&"£>>- 
Wi  =  <3f  =  -|£/^1„, 

The  global  boundary  conditions  at  all  intermediate  nodes  are  zeros  because  there  are  no 
externally  applied  shear  forces  and  bending  moments.  Hence, 


qM  +  Q<a>  =  0  qSi}  +  q!2)  =  0 

•  •  • 

QiN-1]  +  q[n)  =  0  q[n~1]  +  q[n)  =  0 


4It  depends  on  boundary  condition 
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Figure  21:  Adaptive  observer  for  hinged-hinged  beam:  Decay  of  state  error 
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Figure  22:  Decay  of  deflections  under  state  (Riccati)  feedback,  R=0.01 


Figure  23:  Decay  of  deflections  under  state  (Riccati)  feedback,  R=1.64 
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Figure  24:  Decay  of  deflections  under  state  (Riccati)  feedback,  R=100 


Figure  25:  Stabilization  of  a  clamped-free  beam  with  state  (Riccati)  feedback,  R=100 
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:-;ure  26:  Decay  of  deflections  under  output  (k=0.025)  and  state  feedback  (R=100 


Figure  27:  Stabilization  of  a  hinged-hinged  beam  with  state  (Riccati)  feedback 
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imaginary  part 


closed  loop  eigenvalues  under  state  LQR  feedback 


real  part 

Figure  28:  Closed-loop  eigenvalues  for  clamped-free  beam  using  state  LQR  feedback,  8 
elements. 
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Figure  29:  Closed-loop  eigenvalues  for  clamped-free  beam  using  state  LQR  feedback,  16 


elements. 


real  part 

Figure  30:  Closed-loop  eigenvalues  for  clamped-free  beam  using  state  LQR  feedback,  32 


elements. 
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Figure  31:  Closed-loop  eigenvalues  for  clamped-free  beam  using  state  LQR  feedback 
elements. 
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Figure  32:  Maximal  real  part  of  closed-loop  eigenvalues  with  LQR  feedback. 
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Figure  33:  Maximal  real  part  of  closed-loop  eigenvalues  with  LQR  feedback. 
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Decay  of  deflections  under  output  feedback,  kf2=0.02 
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Figure  34:  Decay  of  deflections  under  output  feedback  with  feedback  gain  k=0.02,  the 
5th  modal  shape  function. 
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Figure  35:  Decay  of  deflections  under  output  feedback  with  feedback  gain  k=0.025,  the 
4th  modal  shape  function. 
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Figure  36:  Decay  of  deflections  under  output  feedback  with  feedback  gain  k=0.025,  the 
5th  modal  shape  function. 
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Figure  37:  Decay  of  deflections  under  output  feedback  with  feedback  gain  k=0.03,  the 
4th  modal  shape  function. 
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Decay  of  deflections  under  output  feedback,  kf2-0.05 


Figure  38:  Decay  of  deflections  under  output  feedback  with  feedback  gain  k=0.05,  the 
5th  modal  shape  function. 
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Figure  39:  Decay  of  deflections  under  state  feedback,  R=100,  the  4th  modal  shape  func¬ 
tion 
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Figure  40:  Closed-loop  eigenvalues  of  hinged-hinged  beam  with  state  LQR  feedback, 
R=100,  8  elements. 
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Figure  41:  Closed-loop  eigenvalues  of  hinged-hinged  beam  with  state  LQR  feedback, 
R=10,  8  elements. 
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Figure  42:  Closed-loop  eigenvalues  of  hinged-hinged  beam  with  state  LQR  feedback, 
R=5,  8  elements. 
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Figure  43:  Closed-loop  eigenvalues  of  hinged- hinged  beam  with  state  LQR  feedback 
R=5,  16  elements. 
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real  part 

Figure  44:  Closed-loop  eigenvalues  of  hinged-hinged  beam  with  state  LQR  feedback, 
R=5,  64  elements. 
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